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SUMMARY 


The application of a new, general, potential flow computational technique to the 
solution of the subsonic, three-dimensional flow over wings with leading-edge vortex 
separation is presented. The method is capable of predicting forces, moments, and 
detailed surface pressures on thin, sharp-edged wings of rather arbitrary planform. (It 
should be noted that operational limitations related to doublet panel spacing and 
density requirements, behavior at large planform breaks, convergence characteristics, 
etc., have yet to be extensively explored.) The wing geometry is arbitrary in the sense 
that leading and trailing edges may be curved or kinked, and the wing may have 
arbitrary camber and twist. The method employs an inviscid flow model in which the 
wing and the rolled-up vortex sheets are represented by piecewise continuous, 
quadratic doublet sheet distributions. The Kutta condition is imposed along all wing 
edges. Strengths of the doublet distributions, as well as shape and position of the free 
vortex sheet spirals, are computed in iterative fashion, starting with an assumed initial 
sheet geometry. 

The numerical method has been written in FORTRAN 2.3 language for the CDC 
6400/6600 digital computers. To demonstrate the validity of the method, selected planar 
wing geometries have been analyzed in incompressible flow over a wide range of angle 
of attack. These geometries include: (1) delta wings of differing aspect ratios, (2) an 
arrow wing, and (3) a gothic wing with swept trailing edge. Wing normal force 
coefficients, load distributions, and detailed surface pressure distribution are presented 
and compared, whenever possible, with theoretical results of the leading-edge suction 
analogy and with experimental data. Good agreement with the available theoretical and 
experimental data is demonstrated. 

During the contract period only a limited number of check cases were executed. It has 
since been found that more general wing geometries, particularly ones with kinked 
leading edges, might experience convergence problems. Problems of this type can be 
caused by inadequate panel density, a poor initial guess for the fed sheet geometry, and 
inappropriate modeling at planform breaks where a second free vortex is sometimes 
present. Ongoing work is aimed at improving and clarifying the operational limits of 
the program. 



INTRODUCTION 


The flow at the leading and tip edges of a swept wing with sharp edges separates at 
moderate to high angles of attack, producing vortex sheets that roll up into strong 
vortices above the upper surface of the wing. The formation of these vortices is 
responsible for the well-known nonlinear aerodynamic characteristics exhibited over the 
angle-of-attack range (fig. 1). 

A description of the physical aspects of the flow phenomenon is provided in references 1 
and 2. A line of attachment (point A in fig. 2) is found inboard of the leading and tip 
edges on the lower, or windward, side of the wing. The crossflow velocity components 
inboard of this line are small, so the flow is swept rapidly downstream in this region. 
The flow velocity outboard of this line of attachment is of the order of the freestream 
speed, and streams outward from the leading edge. A similar outward flow occurs on the 
upper surface of the wing with about the same velocity but with a different direction. A 
shear layer is convected from the leading and tip edges to a position above the wing 
under the influence of its own vorticity, where it rolls up to form the principal vortex. 
The air outboard of the shear layer flows first upward, then inboard and around the 
principal vortex, and finally downward toward the wing, forming a line of reattachment 
on the upper surface (point B in fig. 2). The crossflow velocity inboard of this line of 
attachment is small, and hence the air in this region flows rapidly downstream. 

Just outboard of the upper surface attachment, the flow is directed toward the leading 
and tip edges. It accelerates before passing under the principal vortex, and then slows 
as it approaches the edge. The boundary layer sometimes separates because of the 
adverse pressure gradient. The flow immediately reattaches, forming a secondary vortex 
in the opposite sense to the principal vortex, and the flow then continues outward to 
form the shear layer off the leading edge. A weaker tertiary vortex has also been 
observed in the boundary layer near the edge. The secondary and tertiary vortices, for 
simplicity, are not shown in figure 2. Experimental studies of the principal shear layer 
indicate that its form is relatively independent of Reynolds number. 

The rolled-up part of the vortex sheet consists of three regions: (1) an outer region in 
which the distance between turns is large compared to the diffusion distance \fvt, 
where v is the kinematic viscosity and t is the time in which the fluid element in the 
sheet has traveled from the wing leading edge; (2) an inner region where the distance d 
between turns is of the same order of magnitude as the diffusion distance; and (3) an 
inner viscous core that is very small and represents only about 5% of the diameter of 
the vortex sheet. In the outer region, where d > y/v t, convection is the dominating 
influence. In the second region where d is of the order of y/vt, the vortex layers are less 
defined but convection is still important, while diffusion dominates the central 
viscous core. 

The leading-edge suction analogy described in references, 3, 4, and 5 provides a method 
suitable for calculating the magnitude of the nonlinear vortex lift on a rather broad 
class of wing planforms. Polhamus (ref. 3) reasoned that the normal force needed for the 
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Figure 1 .-Leading-Edge Vortex Flow 
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Figure 2.— Cross Section of Flow Over Delta Wing With 
Leading- Edge Vortex (Crossflow Mode!) 
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flow around a leading edge to reattach to the wing is equivalent to the leading-edge 
suction force necessary to force the flow to be attached to the leading edge in an 
unseparated condition. The unseparated leading-edge suction force is calculated, and 
then rotated normal to the wing to obtain the lift contribution of the leading-edge 
vortex. The total wing lift computed by this method agrees well with experimental data, 
but the leading-edge suction analogy does not give flow-field details or detailed surface 
pressure distributions. 

Several attempts have been made in the past to theoretically predict detailed pressure 
distributions and flow fields about swept wings with leading-edge vortex separation. 
Most of these past methods were limited to slender configurations in which considerable 
simplification occurred because the problem could be reduced to a solution of Laplace’s 
equation in the crossflow plane, for which conformal mapping became a powerful tool. 
Smith (ref. 6) developed the best known method of this type by improving the work done 
earlier in collaboration with Mangier (ref. 7). Assuming conical flow, which is 
approximately valid near the apex of the wing, he was able to predict qualitatively the 
type of pressure distributions that had been observed experimentally. Those pressure 
distributions (fig. 3) exhibited a vortex-induced pressure peak at about 70% local 
semispan of the wing. Toward the trailing edge, Smith’s method overpredicted the 
experimental load distribution by a considerable amount, the reason being that his 
conical theory did not satisfy the Kutta condition at the trailing edge. Figure 3 shows 
qualitatively such a comparison of Smith’s theory with experiments and also, for 
illustrative purposes, spanwise pressure distributions from linear lifting surface theory 
(ref. 8) and from R. T. Jones’ slender wing theory (ref. 9) at two chordwise stations of a 
delta wing. This figure (supplied by Mr. Blair Gloss of NASA Langley), shows clearly 
that none of these theories can even approximately predict aerodynamic load 
distributions of wings with leading-edge vortex separation, and demonstrates the need 
for an accurate prediction method for this type of flow phenomenon. 

Several authors have published prediction methods based on the vortex lattice 
technique, among them Rehbach (ref. 10), Mook and Maddox (ref. 11), and Kandil, 
Mook, and Nayfeh (refs. 12 and 13). Sufficient data on pressure distributions have not 
been available from these references for a thorough comparison with the method of this 
document. 

The method described in this report models the flow over wings with leading-edge 
vortex separation as an inviscid and irrotational problem. The method is completely 
three-dimensional and capable of predicting detailed pressure distributions as well as 
overall wing-load coefficients. The wing, wake, and primary vortex sheet are modeled as 
thin surfaces. No attempt is made to model secondary flow separation or multiple 
primary vortex sheets. Special attention is paid to the modeling of the viscous core of 
the rolled-up vortex spiral within the framework of inviscid flow theory. This viscous 
core model is termed the "fed sheet model” in this report. A very simple treatment of 
the core region is adopted, justified mainly because the diameter of the viscous core 
region is small compared to its distance from the wing. The fed sheet model described in 
this report, although adequate for most wing geometries and flight conditions, can be 
improved and must, therefore, be viewed as an interim solution. 


5 





Figure 3.— Load Distribution on Delta Wing Given by Earlier Theories 
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An advanced aerodynamic panel method (ref. 14) is used to solve the mathematical 
model problem describing the leading-edge vortex flow. This new method uses source 
and doublet panel networks as basic building blocks from which a great variety of 
boundary value problem representations can be composed. In the present work, only the 
doublet networks are used to simulate the thin surfaces under consideration, but the 
computational algorithms (e.g., source networks) are available for extending the method 
to handle wing thickness and fuselage simulation in the presence of free vortex flows. 

A short summary of the method described in this document has previously been 
published (ref. 15). 
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SYMBOLS 


[A] matrix of aerodynamic influence coefficients 

Ajj vector in the near plane of a panel surface 

Afcj aerodynamic influence coefficient 

AR aspect ratio 

a parameter of Smith’s conical solution 

[B] matrix, defined by equation (42) 

Bjj vector in the near plane of a panel surface 

b wing span 

[C] matrix, defined by equation (23) 

[Co] matrix, defined by equation (41) 

Coo, Co£, Co-ij, coefficients, defined by equation (60) 

Coff Co**, 

CoT)T) 


cl 


c m 


cn 


Cp 


Acp 


c r 

c ref 

D ( i ) 

[D] 

[DK] 

E 


lift coefficient 

pitching moment coefficient 

normal force coefficient 

pressure coefficient 

jump in pressure coefficient 

root chord 

reference length 

correction term of Jacobian 

matrix, defined by equation (19) 

matrix, defined by equation (40) 

function, see equation (46) 

unit vectors of local panel coordinate system 
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F 

— ► 

F 

f<i) 

G 

J 

[K] 

K(P,Q) 

li,m 

M 

M P 


N 

N d 

N f 

N fs 


N w 


P(f,H,t) 

R 

r*PQ 


function, see equation (46) 
force 

vector, defined by equation (52) 
function, see equation (46) 

Jacobian 

matrix, defined by equation (37) 
potential of an elementary doublet 

chord length of panel segment in transverse geometry cut 
number of panels on one-half of configuration 
pitching moment 
freestream Mach number 

number of doublet parameters in neighborhood of panel 
number of doublet parameters on one-half of the configuration 
normal force 

number of free-sheet panels on one-half of the configuration 
one-half of the number of wing panels 
panel normal vector 
panel corner point 

vector defining panel corner-point location 

panel center 

position of field point 

position of elementary doublet 

residual 

vector from elementary doublet to field point 
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S\v 
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Uk 
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V 


n 


V s 


w k 

x,y,z 


x 0 , yo> zo 


Xp 

Ax 

a 

r 

y 

V 

7 

A 

5 (i) 

«U 

0 


panel area 
wing area 

local wing semispan 
time 

component of freestream velocity 

freestream velocity 

velocity 

normal velocity component 
average sheet velocity 

velocity induced by a panel doublet distribution 
weights 

wing fixed Cartesian coordinate system 
coordinates of panel center 
pitch axis 

distance between transverse geometry cuts in x-direction 
angle of attack 

strength of line vortex along terminated edge 
semiapex angle of delta wing 
constant in Quasi-Newton method 
vorticity 

jump or perturbation of variable 
scaling parameter 
Kronecker delta 

panel inclination in transverse cut 
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A 

sweep angle 


strength of panel doublet distribution 

Me 

doublet parameters along edges of networks 

Mj 

doublet parameters on one-half of configuration 

Mk 

doublet parameters in neighborhood of panel 

Mr 

doublet parameters not located along edges of networks 

M0> M£> Mt b 

coefficients of panel doublet distribution 



V 

kinematic viscosity 

<b 

velocity potential 

{,v,t 

local panel coordinates 

Subscripts 


LE 

leading edge 

L 

left side of configuration 

1 

lower side of sheet 

R 

right side of configuration 

TE 

trailing edge 

u 

upper side of sheet 

x,y,z 

components in wing fixed Cartesian coordinates 

t,v,l 

components in local panel coordinates 

OO 

freestream value 

Superscripts 


D 

difference of values across the sheet 

(i) 

iteration number 

S 

average values of both sides of the sheet 

T 

transpose of a vector 
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THEORETICAL METHOD 


BOUNDARY VALUE PROBLEM 

The flow over wings with leading-edge vortex separation is modeled on the following 
assumptions: 

• The wing is a thin surface (camber surface) in steady, symmetric flight. 

• The flight Mach number is low subsonic; i.e., the flow can be regarded as 
incompressible. 

• Vortex lift is produced by a symmetric pair of primary vortex sheets separating 
from the leading edges of the wing. 

• Secondary flow separation and boundary layer displacement effects are neglected. 

• Multiple primary vortex sheets do not form. 

• The flow surrounding wing and vortex sheets can be regarded as inviscid and 
irrotational, i.e., as potential flow. 

• Vortex breakdown does not take place. 

• The diameter of the viscous core of the rolled-up primary vortex sheet is small by 
comparison with its distance from the wing surface and by comparison with a 
characteristic lateral dimension of the wing. 

The essential elements of this inviscid, irrotational, and incompressible flow model are 
the wing, wake, primary vortex sheet (termed free sheet), and the so-called fed sheet 
representing the core of the rolled-up vortex sheet (fig. 4). 

The flow is governed by Laplace’s equation for the velocity potential, 

v 2 <b = 0 ( 1 ) 

The velocity, V, in turn is obtained from 

V = V4> (2) 

The boundary conditions imposed on the flow model are: 

• The flow must be everywhere parallel to the wing surface, i.e., 

n • V S = 0 (3) 
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The symbol n stands for the unit normal vector of the considered surface. The 
superscript S of the velocity indicates an average sheet velocity, defined by 

V S = '/ 2 (V U + vj) (4) 

— ► — ► 

where V u > Vi are the velocities on the upper and lower sides of the sheet surface, 
respectively. 

• The free sheet and the wake cannot support a pressure differential. This statement 
expressed in terms of the jump in pressure coefficient across the sheet reads 

Ac p = 0 (5) 

In addition, free sheet and wake are stream surfaces; i.e., they must satisfy 

n-V S = 0 (6) 

• Kutta conditions are imposed along the leading, side, and trailing edges of the 
wing in the presence of free sheets and wake emanating from these edges. 

The Kutta condition need not be treated as a separate boundary condition. Satisfaction 
of the stream-surface boundary conditions (3) and (6) of wing, free sheet, and wake will 
result in a smooth flow off the wing edges. The additional requirement of zero pressure 
jump across the free sheet and wake in the immediate vicinity of the wing edges is 
already stated by equation (5). 

• The freestream is undisturbed at infinity. 

• The fed sheet (fig. 4) is an entirely kinematic extension of the free sheet. Its 
dimensions are taken from a conical flow solution used as the "initial guess” in the 
present iterative solution procedure applied to the fed sheet. This is a simplified 
model of the true physical vortex core region, which is viscosity dominated. The 
tacit assumption in this model is that the boundary conditions applied to the free 
sheet are sufficient to adequately position the fed sheet. Appendix A explains why 
this model is used instead of one in which the fed sheet is required to be force-free. 

GEOMETRY DEFINITION 

GENERAL DESCRIPTION 

The geometry of the configuration (fig. 5) consists of four elements: wing, free sheet, 
wake, and fed sheet. 

The wing is represented by the mean geometric surface (camber surface) with an 
arbitrary distribution of camber and twist. The wing must have a pointed upstream 
apex. Leading and trailing edges of the wing may be curved or straight, and the latter 
may be swept forward or backward. The panel scheme described replaces curved wing 
edges with straight-line segments. 
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The free sheet, a thin surface attached to the leading edge of the wing, extends in the 
streamwise direction from the apex of the wing to the point of maximum local wing 
span. The shape of the free sheet in a transverse cut through the configuration 
approximates a spiral. Approximately one-half of the first turn of the spiral provides a 
sufficient model in most cases. 

The wake is a thin surface attached to the trailing edge of the wing downstream of the 
point of maximum span and includes the extension of the free sheet in the downstream 
direction (fig. 5). The wake geometry is generated by straight lines parallel to the x-axis 
of the wing-fixed Cartesian coordinate system. 

The trace of the fed sheet surface in a transverse cut (y,z-plane) through the 
configuration is a straight line connected to the last point of the free sheet and is 
approximately perpendicular to the free sheet. 

The geometry of the configuration is symmetric with respect to the x,z-plane. 

PANEL SCHEME 

The geometry of the configuration is defined by panel corner points. A wing-fixed 
Cartesian coordinate system is employed to specify corner-point positions. The origin of 
this x,y,z-coordinate system (fig. 5) coincides with the apex of the wing. The x-axis 
(positive in downstream direction) is parallel to the wing root chord. 

The geometry of the configuration is paneled as follows. The x-coordinates of transverse 
planes cutting the configuration (fig. 6) are specified. The distance Ax between cuts need 
not be equal. These transverse cuts define the x-coordinates of all panel corner points. 
The wake region is not cut, but is bounded by a transverse cut far downstream of the 
wing. The computer code places this last cut at 50 wing root chords downstream of 
the wing. 

Corner points of the panels are then defined along the geometry contours of each 
transverse cut. This must be done such that the number of wing panels is the same in 
each spanwise row of panels. The same rule must be observed when paneling the free 
sheet. Figure 7 shows the panel corner points in a typical transverse cut through wing 
and free sheet. The fed sheet is represented by a single panel in the y,z-plane. 

The paneling of the wake is completely determined by that of the wing, free sheet, and 
fed sheet. A single wake panel is used in the downstream direction. The wake panels 
are generated by straight lines parallel to the x-axis. 

Examples of the wing panel scheme are shown in figure 8. Because of symmetry, only 
one-half of the configuration needs to be paneled. Curved wing planforms are 
approximated by straight-line segments. (The algorithms extending the numerical 
representation to encompass curved panels are now available, ref. 14, but were not 
implemented in the present work.) Figure 9 displays a spanwise row of quadrilateral 
free-sheet panels between two adjacent transverse cuts. An example of the wake 
paneling is given in figure 10 in which the wake panels downstream of the free sheet 
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Far downstream 


Figure 6. — Transverse Cutting Planes of Configuration 
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Figure 7.— Panel Corner Points in Transverse Cut of Configuration 
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Free-sheet panel 




Figure 10-Example of Wake Paneling (Not Complete) 
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are not visible. It can be noted in figure 10 that in the spanwise direction the number of 
wake panels is not always the same as the spanwise number of wing panels. 

PANEL GEOMETRY 


Consider the single panel shown in figure 11 located between two transverse cuts, Xj and 
xj +1 . The same type of panel geometry is used for the paneling of the wing, free sheet, 
fed sheet, and wake. Triangular panels such as those at the apex of the wing are 
special cases of the shown quadrilateral panel shape in which two corners coalesce. 


The four corner points of a panel, Pj j; Pjj+x; Pj+jj; Pi+ij+i, are used to define a 
planar panel surface that approximates the nonplanar geometry surface. This planar 
surface is the so-called near plane defined by and containing the vectors Ajj and Bjj (see 
fig. 11). The latter are written in terms of vectors that define the panel corner point 


locations as 


'iJ 


P ij + P iJ+l 


- p i+lJ- p i+lJ+l 


( 7 ) 


B ij P ij + P i+l,j P i,j+1 P i+l,j+l 

The panel normal vector n at the center Pq^ of the (ij)th panel is then 


n = 




n = (n x , n y , n z ) 


( 8 ) 


Projecting the four corner points onto the near plane and connecting the new points by 
straight lines defines the panel geometry. All original corner points turn out to be 
equidistant from their projections onto the near plane. 


In addition to the global Cartesian x,y,z-coordinate system, a local panel coordinate 
system (fig. 11) is introduced whose ^-coordinates are defined as follows. The 
^-coordinate points in the direction of the panel normal; the 7j-axis is perpendicular to £ 
and to the x-coordinate of the global coordinate system. The £-axis (positive in 
downstream direction) forms, with 17 and £, a right-handed Cartesian coordinate system. 
The and x-axes are in general not parallel. The unit vectors of the local panel 
coordinate system are denoted by e£, e^, e£. 


The transformation from x,y,z- to £, ^-coordinates is performed by 


r ^ 

% 


n y 2 + n z 2 

n x n y 

n x n z 


r 

X 

- x 0 


■yV + n z 2 

\P + n z 2 


* V 

► - 

0 

n z 

.... n y 

i 

y 

-Vo 


y n y 2 + n z 2 1 

/ n y 2 + n z 2 


r 

j 


" x 

n y 

n z 


z 

_z 0. 


where the coordinates of the panel center Pq..= ( x o>yo> z o)ij are average values of the 
original panel corner-point coordinates. 
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Note: All original corner points turn out to be equidistant 
from their projections onto the near plane 


Figure 1 1 .-Single-Panel Geometry 
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NUMERICAL METHOD 


INTRODUCTORY REMARKS 

The numerical solution of the boundary value problem is accomplished by a panel-type 
influence coefficient, method. The numerical method formulated in this section is 
entirely based on an advanced panel method for the solution of Laplace’s equation 
subject to Neumann, Dirichlet, or mixed boundary conditions (ref. 14). Only certain 
features of this new method that center about the use of doublet distributions to 
represent the various elements of the adopted aerodynamic model are employed. 

The surface of the configuration is divided into panels, as explained in the section on 
"Geometry Definition.” Doublet distributions dependent on a finite number of unknown 
doublet parameters are defined for each panel. A finite number of control points on the 
surface are selected at which the boundary conditions are satisfied. This discretization 
of the problem leads to an aerodynamic influence coefficient formulation wherein the 
solution is achieved by assembling logically independent networks of doublet panels. In 
this context, a network is defined as a portion of the boundary surface on which a 
certain distribution of doublet strength is specified, together with properly posed 
analysis (Neumann) or design (Dirichlet) boundary conditions. It is always composed of 
a rectangular array of panels that may be input independently of all other networks. A 
network is logically consistent in that it contributes as many equations to the overall 
problem as it contributes unknowns. 

The essential features of the numerical method are: 

• Discrete values of doublet strength are assigned to certain standard points on each 
network. A local distribution of surface singularity strength is obtained by fitting a 
quadratic doublet form to these discrete values in an immediate neighborhood by 
the method of least squares. 

• Certain standard points on each network are assigned as control points. These 
points include panel center points as well as edge abutment downwash points. The 
latter serve to automatically impose standard aerodynamic edge conditions, e.g., 
the Kutta condition, zero potential jump at thin edges, continuity of singularity 
strength across abutting networks, etc., producing logical independence for each 
network. In all cases, the number of boundary conditions on each network coincides 
with the number of assigned surface singularity parameters. 

• Two expansions are employed in the calculation of the aerodynamic influence 
coefficients: a near-field expansion and a far-field expansion. All resultant 
integrals are evaluated in closed form using recursion relations, which contain the 
fundamental logarithm and arc tangent transcendental terms that appear in 
constant-strength panel techniques. 
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DOUBLET PANEL 

The distribution of doublet strength ^(£, 17 ) on a panel is quadratic; i.e., 

= Mq + t? + y £ 2 + ^ -^t? 2 ( 10) 

The six coefficients, /xq, /x^, /x^, /x^ are not assumed independent; rather they 

are a linear combination of an independent set of doublet parameters /xj. The 
parameters /xj are the singularity strengths at a set of discrete points on the network 
surface. The linear relationships between the six coefficients and independent doublet 
parameters are later determined by the method of weighted least squares. 

Consider now the velocity induced by a distribution of elementary doublets over a 
single panel S at a field point P = (£, v>0- 

v(P )=ff m(Q) V p K(P,Q)dS 
S 

The integration is to be carried out over the panel area S with Q 
position of the elementary doublet on the planar panel. Further, 

d S = d d 7?i i\2) 


(ID 

= (£1 j 1 ? 1 ) denoting the 


K(P,Q) =47 


1 r PQ • n Q 
r PQ 3 


(13) 


IfpQl = r pQ 

where K(P,Q) is the potential of an elementary doublet. The vector ?pq is 


rpQ = (£-Si)e£+0? + 


(14) 


and points from the position Q of the elementary doublet to the field point P. The vector 
Sq is the panel normal i?q = (0,0,1) for planar panels. Applying the gradient operator 


3 — ► , 3 — ► , 3 — ► 
V P = d £ 3rj e r? df e f 


(15) 


to K(P,Q) yields 


V P K(P,Q) 



3? PQ 

r PQ 5 


(? PQ ‘ 


(16) 


Appendix B contains details of the integration of equation ( 11 ). The analysis 
distinguishes between: 


1 . 


Field points P that are located in the vicinity of the panel S, in which case no 
approximation to the integrand for flat panels is used. Then 


1 



U 


K 



K = 3,5 (17) 
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2. Field points located a sufficient distance away from the panel surface to justify the 
following approximation to rpQ~ K 


— 4 = U K + KU K+ “(££ j + 7717 ] ) ; K = 3,5 
r PQ 


U 


1 

Vs 2 + *? 2 + r 2 


(18) 


This approximation is applied when U is greater than 2.5 times the maximum 
panel diagonal. 


The result of the integration of equation (11) can be written in matrix form as 


f a 

V; 


IT 


= m 


V 0 

T 

v. y 


(19) 


where v*, v^, are the components of the perturbation velocity v in local panel 
coordinates of the inducing panel. 


NETWORKS 

The configuration is modeled aerodynamically by five networks termed: 
• Analysis network (wing) 


• Design network (free sheet) 

• Fed-sheet network 


• Wake network 

• Network number 5 

The analysis- and design-type networks are the two basic networks, whereas the 
remaining three are specializations of the analysis network. The last network is needed 
for that wake panel downstream of the last fed-sheet panel. 

Figures 12 and 13 show schematically the arrangements of doublet parameters /Aj and 
the control point locations for the analysis and design networks, respectively. These 
choices of the locations of the doublet parameters ensure a stable computation. For each 
of these networks the number of control points at which the boundary conditions or edge 
conditions are satisfied equals the number of unknown doublet parameters. The edge 
control points serve in the matching of the various networks across their common 
boundaries. These points are slightly removed from the edges for this purpose (10 5 of 
the maximum panel diagonal). 
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An individual panel of both the analysis and the design network is assumed to have a 
quadratic distribution of doublet strength, which is defined by equation (10) in terms of 
six coefficients. A linear relationship between these coefficients and those doublet 
parameters p* on the panel itself and on its immediate neighbors is determined by the 
method of least squares. For each panel S the residual R is minimized. 



W k[M(^k>-M k ] 2 


(20) 


The weights are chosen large (10 5 ) when (£k,i7k) actually lies on the panel S, and 
are unity if not. The number N of the doublet parameters fi k depends on the position of 
the panel within a network. In the case of a panel of an analysis network, N equals 9 
unless the quadrilateral panel and its neighbors are triangular. For most panels that 
are part of a design network, the summation is carried out over 16 doublet points. Only 
panels located along the edges of a design network have fewer neighboring doublet 
points, in which cases N reduces to 9 or 12. 

With given by equation (10), the residual R reads 

1 ^ j _ i 2 

^ = ~2 ^ ^k(^0 + ^|^k + ^r/^k + ~ ^|^k + ^^k^k + "2 - ^k) > 

k=l v 


which is minimized by 


|B-=o #^. = 0 = o 

8 <*0 iv -r, ^iri %t| 


This procedure gives the following intermediate result 
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with 
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The desired linear relationship between the six coefficients of the quadratic doublet 
distribution and the N doublet parameters can then be written as 

R1 


< ► 

* nn j 



(23) 
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where the coefficients of the matrix [C] depend only on the chosen weights W k and on 
the geometry of the panel scheme. 

Figure 14 shows the wake network unwrapped. This network is of the analysis type but 
is simplified assuming that the doublet distribution varies only in the lateral or 
77 -direction. 

Figure 15 illustrates the chosen positions of doublet parameters and control points for 
the fed-sheet network, which is like the wake network of the analysis type. The doublet 
strength is constant in the redirection, and varies quadratically in the direction parallel 
to the terminated edge. 

The number 5 network needed for the wake panel located downstream of the fed sheet 
has a constant distribution of doublet strength. This is a consequence of the choice for 
the doublet distributions of the neighboring fed sheet and wake panels. 

MATCHING OF NETWORKS 

The relative positions of the five networks described in the previous section are 
illustrated in figure 16, which also shows the location of most of the control points, 
particularly those along the edges of networks. The reader will notice that at some 
boundaries edge control points appear in opposing pairs. This is the case along the 
centerline of the wing and along leading and trailing edges of the wing. At other 
boundaries such as the one between free sheet and fed sheet and between free sheet and 
wake, there is only a single point controlling the edge condition. 

— ►o 

The downwash V n = n * V is specified at edge control, points. In order to explain the 
implications of this edge condition, the downwash in the vicinity of the common edge of 
two smoothly adjoining networks is expanded as 


V n (s) = 




In s + regular terms 


(24) 


Here, s is a coordinate tangential to the network surface and perpendicular to the edge, 
A^ e is the jump in doublet strength across the edge, and A(dfi/ds) e is the jump in the 
derivative of doublet strength across the edge. A control point placed near the edge 
(e.g., at s = 0+) requiring that downwash be finite will tend to make A/i. e vanish, i.e. /t e 
continuous across the edge. This can easily be shown by taking the limit 

lim Aju e (s) = lim j sV n (s) - a (^~) slns-s(. ..)[ = 0 (25) 

s-0+ s-0+/ V /e ’ 


Placing a second control point on the opposing panel of an adjoining network will, in 
addition, force A(d/i./ds) e to vanish, thereby establishing continuity of (dfj./ds) e . 


28 



Wing 




Free sheet 



Doublet parameter 
Control point 


Unwrapped 




Wake panel 


Figure 14.— Wake Network (Simplified Analysis Network) 
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Figure 15.— Fed Sheet Network (Simplified Analysis Network) 
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Hence, across the centerline of the wing and across leading and trailing edges, /z as well 
as d/z/ds are continuous, since control points are grouped in pairs along these edges. In 
addition, the derivative of /x must vanish along the centerline to produce a symmetric 
distribution of /z. As a consequence of the particular construction of the wake network, 
the derivative of /x in the freestream direction vanishes along the trailing edge. 

Doublet strength, but not its derivative, is continuous across the common boundaries of 
free sheet and wake, free sheet and fed sheet. Thus it will be seen that the network 
edge downwash points are, in essence, a logical algorithm designed to produce the 
desired degree of continuity in doublet strength across network edges, and are 
functionally quite different from the downwash points appearing in the centers of 
panels. 

AERODYNAMIC INFLUENCE COEFFICIENTS 


In subscripted notation, aerodynamic influence coefficients A^j are defined by 

k= s>i?,-r 

= ( A kj + U k> e k (26) 

j =1,2 ,...N d 


C 

where V is the average velocity at a particular control point, /Zj are the doublet 
parameters of one-half of the configuration, N^ is the number of these doublet 
parameters, and are the components of the freestream velocity in the direction of the 
local panel coordinates 

iL = (U^U^Uj) (27) 

The vectors e^, e^, are the unit vectors of the local coordinate system of that panel on 
which the control point lies. For future reference, equation (26) is rewritten in the form 
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for the components of V in local panel coordinates. 


(28) 


VS = (v^S V/, V r S ) (29) 

— ► c 

In particular, Vj is the component of the velocity normal to the panel surface at a 
control point, i.e., the downwash 


V f S = V n = n • V s 


(30) 
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To derive expressions for the aerodynamic influence coefficients, the analysis proceeds 
by first computing average perturbation velocities at a control point due to a panel 
distribution of doublet strength. The latter is represented by the six coefficients of the 
quadratic doublet distribution (see eq. 10). So the objective is to find equations of the 
form 





(31) 


The perturbation velocity components \£, vfy v^ induced by a panel doublet distribution 
are already known from equation (19). These are also the average perturbation velocity 
components of both panel sides unless the control point is located on the influencing 
panel itself. In this case the average perturbation velocities change to 


r 

, s 


i 



v r 2 y v 

v ^ 
V 

► = < 

vin 

„ s 



v J 


v ? 

V* / 


where y^, y^ are the components of the sheet vorticity y in local panel coordinates. 
Vorticity y* and doublet strength /z are equivalent concepts in potential flow theory 
related by 


7 = nx Vju 


(33) 


Recalling the definition of the panel normal vector n = (0,0,1), the components of y can 
be written in terms of the components of V/z as follows 


d/z _ dfJi 

dr] 



(34) 
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or upon introducing //(£,??) from equation (10) 


Hence, 


1 

2 7 h 


r M 0 > 



1 



2 t £ 

► = [K] < 

Hi 
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^ J 


' J 


( 35 ) 


(36) 


with 


[K] 




(37). 


Combining equations (19), (32), and (36), the average perturbation velocities induced by 
a single panel take the form 



r° l 


[[D] + [SjjKl ) 




(38) 
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The Kronecker delta 8{j indicates that the coefficients of the matrix [SjjK] are zero if 
the control point i does not lie on the influencing panel j(i*j), and that they are given 
by equation (37) if i = j. 


At this point of the derivation of aerodynamic influence coefficients, advantage is taken 
of the symmetry properties of the flow field. The assumption is made that the vortex 
separation is completely symmetric with respect to the x,z-plane of the wing. 
Consequently, the distribution of doublet strength is also symmetric; i.e., each doublet 
panel on the RHS (subscript R) of the configuration has an image of equal strength on 
the LHS (subscript L). Hence, 


However, the velocity perturbation caused by a panel of the RHS at a particular control 
point differs from the perturbation induced at the same point by the image panel. The 
sum of the average perturbation velocities of a panel and its image at a control point on 
the RHS of the configuration takes the form 
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V ? S J 


rn 0 


> = [DK] < 




ri 
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r?r? 


(40) 


with 


[DK] = [D] R + [ 5j jK ] R + [D] L 


Next, the six coefficients of a panel doublet distribution must be related to the doublet 
parameters /xj defined on one-half of the configuration. Recalling that equation (23) is a 
linear relationship of the panel coefficients and the set (i ^ of doublet parameters of the 
immediate neighborhood of the panel, this is easily accomplished by enlarging the 
coefficient matrix [C]. Then 






> = 


= Mb 


v v 


(41) 
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where most of the coefficients of the enlarged matrix [Co]are zero. An indexing system 
explained in the description of the computer code identifies the nonzero coefficients 

°f [C 0 ] ■ 

Introducing equation (41) to (40) yields 


= [B]j/Xjj [B] = [DK] [C 0 ] (42) 


The reader is reminded that equation (42) only gives the velocities induced by a single 
panel and its image, even though {/ij} is the vector of all unknown doublet parameters. 
The perturbation velocities induced by all M panels of one-half of the configuration and 
their images are obtained from 




The matrix [ A] is the matrix of aerodynamic influence coefficients, i.e., 


(43) 


[A] - 


A fj 


A r- 

s j 


which was introduced earlier by equation (28). 

SOLUTION PROCEDURE 

The boundary value problem of wings with leading-edge vortex separation is nonlinear 
because the shape of the separated vortex sheet as well as its strength are unknown. 
The solution procedure must therefore be iterative. Details of the iteration procedure 
are discussed in this section. First, the choice of initial values for the free-sheet 
geometry and the doublet distribution is described. The result is termed "Starting 
Solution” of the iteration. Then a description of the perturbation technique follows and 
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is used to update free-sheet position and doublet strength during each step of the 
iteration, entitled "Update Scheme.” Finally, an outline of modifications of the basic 
update scheme, introduced to ensure stability and economy of the computation and 
known as the "Quasi-Newton Method,” concludes the description of the solution 
procedure. 

STARTING SOLUTION 

An initial guess for the free-sheet geometry is established based on conical solutions of 
Smith (ref. 6). Smith’s results are reproduced in figure 17, which shows the shape of the 
free sheet and the size of the fed sheet for various values of the parameter a. The latter 
is defined by 


a 

tan 7 


(44) 


where a denotes the angle of attack in radians and y is one-half of the apex angle of a 
delta wing. The sheet geometries of figure 17 represent transverse cuts through the 
configuration normal to the wing surface. The y,z-coordinates are nondimensionalized 
by the wing semispan s. The locations of the line vortex along the terminated edge of 
the fed sheet are given for several values of a (0.2 ^ a 3.0) and are connected by a 
dashline. The straight line between the last point on the free sheet and the line vortex 
is the trace of the fed sheet. An example is shown for a = 1.4. 

The shape of each free sheet is subdivided into seven segments of approximately equal 
chord length. The corner points of the segments, as well as the nondimensional 
coordinates of the position of the terminated edge, are contained in the computer code. 

Figure 18 explains how an initial free-sheet geometry is obtained for a nonconical wing 
geometry based on the tabulated data of Smith. For this purpose the assumption is 
made that initially the shape of the free sheet at a particular chordwise station is the 
same as that of a certain delta wing. This delta wing is locally equivalent to the 
considered nonconical wing geometry and is defined as a wing that has the same apex 
position and the same local semispan at that chordwise station where the initial 
free-sheet geometry is to be computed. Thus, the parameter a can be calculated at each 
transverse cut for a given angle of attack and a given angle y=arc tan (s/x). Linear 
interpolation of Smith’s data provides the desired initial free-sheet geometry for a 
chosen number of free-sheet panels. All free-sheet segments of a transverse cut 
(y,z-plane) have approximately the same chord length. 

The described procedure also provides the size of the fed sheet at all geometry defining 
transverse cuts. The lateral length of the fed sheet and its relative position (i.e., local 
angle with the end of the free sheet) to the neighboring free-sheet panel is fixed 
throughout the subsequent iteration procedure; however, it changes absolute position in 
response to movement of the free sheet during the iteration. 

It should be emphasized that Smith’s conical data provide only the initial free-sheet 
geometry. This is a convenient choice and a good guess for wing geometries that are not 
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too different from delta wings. The computed doublet distributions and the sheet 
geometries computed in subsequent cycles of the iteration procedure are, in general, not 
conical. 

An initial guess for the doublet distribution of wing, free sheet/wake, and fed sheet is 
computed by satisfying identically the stream-surface boundary condition 


— ► '“►c 

n • V b = 0 


at all boundary points. This equation is rewritten in terms of aerodynamic influence 
coefficients and doublet parameters as defined in equation (26) 


f y 


+ u> = o 

s j 


where the subscript £ denotes components in the direction of the local panel normal at a 
particular control point i. Introducing matrix notation allows the computation of the 
initial guess of the unknown doublet distribution {jij (1) } from 



(45) 


UPDATE SCHEME 

Beginning with an initial guess for the sheet geometry and the doublet parameters, the 
sheet position and doublet strength must be updated so as to satisfy more closely the 
pertinent boundary condition of the leading-edge vortex problem. For this purpose the 
governing nonlinear boundary condition equations are perturbed, and with the 
assumption of small perturbations, a linear set of equations is derived for the 
perturbation variables. The perturbation technique shall now be discussed in detail. 

The boundary value problem can be written symbolically in terms of the following three 
equations: 


E(/i e ,j Lt r ,0) = 0 

F(p e , M r , d) = 0 (46) 

G(M e ,M r ,0) = 0 

In this notation, the are those doublet strength parameters defined at the edge 
points of networks across which the Kutta condition has to be satisfied. The /z r are all 
remaining unknown doublet parameters. The angle 6 stands for the unknown free-sheet 
position and denotes the panel inclination in a transverse cut or y,z-plane. The above 
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equations state the boundary conditions of the problem and in particular have the 
following meaning: 

E = 0 expresses the stream-surface boundary condition, n * V s = 0, at edge points of 
networks where the Kutta condition is to be satisfied. 

F = 0 symbolizes the boundary conditions Acp = 0 of free sheet/wake and n • V s = 0 of 
the wing. 

G = 0 stands for the stream-surface boundary condition, n * V s = 0, of free sheet 
and wake. 

The above equations are expanded in a Taylor series in which all second-order terms of 
the perturbation variables A//, e , A/x r , A 0 are neglected. 

E (i + ,) = E (0 + ^ iMe+ |^^ + |E A# 

p(i + l) = F 0) + ^ AMe+ |L A , r + |F A9 (47) 

G (i+i ) = G (o + ie. A(ie + 1*. A „ r + |g A <) 


The superscripts (i) and (i + 1) indicate the (i)th and (i+l)th cycle of the iteration, 
respectively. All partial derivatives are known functions of the known values 
of the (i)th iteration cycle. Moreover, the following abbreviations are used 

.E^ i+1) = E( Me ( i+1 ),/x r (i+1) ,0 (i+n ) 
pO+D = F( Me (i+1) ,M r (i+1) ,0 (i+1) ) 

G« . = G(/i e (i) ,M r (i) ,fl (0 ) 


Equations (47) can be written in matrix form as a set of linear equations governing the 
perturbation variables. The assumption is made that all boundary conditions are 
satisfied at the end of each iteration cycle, i.e., 

E (i+1) = 0 

p(i+D = o 

c( i+1 ) = 0 
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In addition, it is assumed that as a necessary requirement for satisfying the Kutta 
condition along wing edges, the stream-surface boundary condition is always satisfied at 
edge boundary points; i.e., E (i) = = 0. This condition is not sufficient for the 

Kutta condition to be satisfied, but by nature of the free-sheet boundary condition, a 
second condition, Acp = 0 just slightly outboard of wing edges, will be satisfied when the 
solution is achieved. Hence, the equations governing the perturbation variables read 
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(48) 


According to the terminology used in Newton-type iteration schemes, the coefficient 
matrix is called the Jacobian. The coefficients of the Jacobian are calculated with the 
following assumptions: 

• Changes. of the aerodynamic influence coefficients due to changes of the panel 
inclination 0 are neglected during each cycle of the iteration. 

• The vector normal to the free sheet at control points in the vicinity of wing edges 
is not affected by changes of 9 . 

• The length of panel segments in transverse geometry cuts, x = constant, does not 
change during the iteration. (This is a fixed constraint throughout the iteration 
that serves to fix the transverse length of the free sheet.) 

• The free-sheet geometry is updated such that panel corner points remain in their 
initial transverse cuts. 


Details of the derivation of the Jacobian are contained in appendices C and D. The 
listed assumptions imply, in particular, that the coefficient 



The number of equations in (48) can be reduced by eliminating A as follows: 



( 49 ) 
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Upon introducing this expression to (48) one obtains 



( 50 ) 


QUASI-NEWTON METHOD 

This scheme (ref. 16) is used to solve for doublet strength and free-sheet position. 


Represent the quantities of equation (50) symbolically as Ax = (A/i r ,A0), the coefficient 
matrix (Jacobian) as J, and the right-hand side as (-0. Equation (50) becomes, 


JAx = -f (51) 

These equations are solved iteratively. Represent the ith iteration by superscript (i). 
The scheme proceeds to find the corrections Ax (i) from the equation 

j(i) Ax^ = -fW (52) 


and forms the new approximate solution (next iterate) 

x (i+l) = xW + fiW AxW (53) 

where J (i) = J(x (i) ), f (i) = f(x (i) ) and 5 (i) is a scaling parameter to limit the step size of 
the correction vector. The Jacobian at x (i+1) is obtained by using the following update 
formula 


j(i + l) = j(i) + D<i) (54) 

where 

D (i) = <f0 +1 > -fCO- j(0 ax^JaxW 7 

Ax (i)T Ax (i ^ 

In this way, there is no need to reevaluate the n 2 partial derivatives for the Jacobian at 
every iteration. The superscript T denotes the transpose of a vector. 

Since the aerodynamic influence coefficients form an essential part of the method, a 
procedure of generating new aerodynamic influence coefficients only after every five 
iterations is included in the iterative scheme. This approach can help to reduce the 
overall computing costs. 
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The scaling parameter 8 (i) is introduced to alleviate the problem of overshoot in the 
classical Newton scheme. For each iteration cycle, the following criteria are used to 
determine 8 (l) 


0 < 5 (i ) < 1 


and 


§(0 HAx^ll < 7 llx^H 

where y is a predetermined quantity {y = 0.1 is presently used in the computer code), 
and || || is the Euclidean norm representing the length of a vector. In addition, a 
halving process of the scaling parameter 8 (i) is applied to ensure the inequality 

||f(i+D|| < ||f(i)|| 

The quality of the solution is monitored by examination of the residual defined by 



(55) 


where k ranges over all appropriate boundary condition points. 

The edge doublets fi e are updated at the end of each iteration cycle. 

WING LOAD AND SURFACE PRESSURES 


Once a converged solution has been calculated the pressure jump across the wing and 
the pressures on upper and lower wing surfaces are obtained from the following 
equations, which are derived in appendix E. 


Ac P = V S • V/z 
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= 1 ^ (v s • V s - V s • Vp + jVp • Vp) 
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(56) 


(57) 


(58) 


The average velocity V s = (\^ S , Vtf, v£ S ) is calculated from equation (28) with the final 
converged distribution of the doublet parameters pj. The gradient of p can be 
determined from 


VP - (P£ + ^ + p^ r} )e^ + + p^| + Vrmrfri (59) 
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which is derived by applying the operator 

V - ttt: Qy + T7 — e„ + zr-r ty 

3£ 5 dr] V s 

to the quadratic distribution of doublet strength yx(£,Tj) of the panel (see eq. 10). The 
coefficients etc., of the doublet distribution, in turn are linearly related to the 

distribution of doublet parameters by equation (41). This equation can be written as 
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(60) 


where Cqq, Cof, Cqtj etc., denote the rows of the coefficient matrix [Cq]. With this 
notation becomes 


- Jcg£ + + 7 ? C 05rj]j^j| e | + [ C 0r? + ^ c 0^7? + T?C 0i7T?])^j I e r? (61) 


Having calculated Acp from equation (56), the normal force coefficient of the wing c^ 
defined by 

Np 

C N = ~ D ^ (62) 

i u ~ 2 % 

becomes 
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C N = Zj Ac P/"z S) j 


(63) 
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The various symbols have the following meaning 
N f normal force 
S\y total wing area 
S panel area 

n z component of normal vector in z-direction 
N w number of wing panels on one-half of the configuration 
The pitching moment coefficient c m , defined by 

Mp 

c = 

m p 9 

*2 U ~“ c ref S W 

takes the form 

2 % r i 

Cm = S Ac Pjh (x o- x p> s Jj 

where 

Mp pitching moment 
xq x-coordinate of panel center 
xp pitch axis 
c re f reference length 


( 64 ) 


( 65 ) 
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COMPUTER PROGRAM USE 


PRACTICAL INSTRUCTIONS 

Some practical hints are given in this section to assist the user in preparing a computer 
run. The code is written in FORTRAN 2.3 for the CDC 6600 digital computer and is also 
operational on the CDC 6400 system. The program occupies 120 000 octal locations of 
central memory in an overlay structure, and uses eight disk files. Many cases have been 
executed to check out the code and to gain confidence in the solution procedure. Most 
computations were highly successful, but there were a few that did not converge to a 
solution within the allotted computer time. These were mostly wing geometries with 
kinked leading edges. The user must be warned that he might occasionally encounter a 
case which, in spite of careful preparation of the input data, does not converge as 
rapidly as desired. If this occurs, numerical experimentation with different paneling 
and/or different choices for the initial geometry might hasten convergence. 

Number of Unknowns : The code is presently limited to 130 unknown doublet 
parameters p j and panel inclinations 6 (in sum total). The doublet parameters ft e 
located along the centerline and the edges of the wing are not counted in the sum. 
There are as many unknown doublet parameters p T on wing and free sheet as there are 
panels, and the number of unknown angles of the free sheet also equals the number of 
free-sheet panels. Hence, the total number of unknowns is equal to the number of wing 
panels N\y plus twice the number of free sheet panels Np<§, and must satisfy 

N w + 2N ps < 130 (66) 

An example of a panel layout is given in figure 19. Table 1 provides additional 
information for this particular case. 

Number of Wing Panels: At least 25 wing panels should be used. Cases with sparser 
wing paneling probably converge to a low-quality solution. The user should also recall 
that the number of panels in all spanwise rows of wing panels must be the same. 

Relative Panel Size: The numerical scheme is not sensitive to strong variations in panel 
size. Convergence difficulties, however, might be encountered when using a small 
free-sheet panel next to a large wing panel along the leading edge. 

Number of Iterations and Updates of Aerodynamic Influence Coefficients: Sixteen 
iteration steps should be sufficient to obtain a converged solution in most cases. The 
program automatically updates the AIC’s after every five iteration steps, so that during 
15 cycles the aerodynamic influence coefficients are calculated four times. 

Monitoring of Convergence: The residual defined by equation (55) is a measure of how 
closely the boundary conditions are satisfied. The residual must decrease rapidly in an 
iteration with good convergence characteristics. In addition, the behavior of the wing 
normal force coefficient cjsj as well as the pressure jump Acp across the free sheet should 
be monitored to judge the quality of the solution procedure. 


46 




47 









Table 1 .-Arrow-Wing Data 


Network 

Panels 

Control points 
at 

panel center, 
doublets 

Edge control 
points 

0 

Wing 

30 

30 

14 

0 

Free sheet 

48 

48 

7 

48 

Fed sheet 

6 

0 

8 

0 

Wake 

11 

0 

13 

0 

Network number 5 

1 

0 

1 

0 

Sum 

96 

78 

43 

48 


Total number of unknowns = 78 + 48 = 126 


( 
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Computer Time: The code uses most of the execution time for the computation of the 
aerodynamic influence coefficients. In order to obtain a first estimate of the total 
solution time, the following procedure is suggested. First, a trial run with one iteration 
cycle and one AIC update should be made. Multiplying the computation time of this 
trial run by the total number of AIC updates and adding approximately 25% to it gives 
an estimate of the total execution time. 


INPUT DATA PREPARATION 

The keywords with the $ sign are important for identifying the corresponding input 
data values. Only the first three characters of the keywords following the $ sign will be 
needed in the program. All data cards with numerical values use format 6E10.0. Table 2 
shows an example for the printed card image of the input data. 




Variable 


Card 

Card 

Name 


Number 

Column 

(Default Value) 

Description/Comment 

i 

1-5 


$CASE 

2 

1-70 


Title information 

3 

1-70 


User identification 

4 

1-16 


$ALPHA (DEGREES) 

5 

1-10 

ALPHAD 

Angle of attack in degrees 

6 

1-13 


$ ASPECT RATIO 

7 

1-10 

AR 

Wing aspect ratio (defined as 4s/x): 
this value will be used for delta 
and arrow-wing preprocessors 

8 

1-16 


$TRANSVERSE CUTS 

9 

1-10 

TRAN 

Number of transverse cuts of 
the wing network (=s 10) 

10 

1-60 

X(I) 

x-coordinates of transverse cuts 
for the wing network 

11 

1-14 


$SPANWISE CUTS 

12 

1-10 

SPAN 

Number of spanwise cuts of the 
wing network (^ 10) 

13 

1-60 

YSP (I) 

Percent values for spanwise cuts 



(0.0) 

(100% = 1); this card is intended 
for use with preprocessor 

14 

1-15 


$CENTERLINE 

15 

1-10 

CTRA 

Number of transverse cuts 
along centerline 
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Table 2. —Card Image of Input Data 


■ A C0MPU T F» PPOGPA** 

JPQP 

4 tmb F e niMFN«*IONU FOLU T TON OF ^LOWS OVFP HI NGF 
WITH LEAOTNG FPGF VCR^EX SEPARATION 



- LIFT OF TNPU T OATA CARPS - 



C A °D IMAGES 




FCAFF 



2 

I3CAS r - P?EnTrTT N G WING Af90 n YN A* t q 

LPACS 

F^p 4 OFL-A HTnO 

3 

IDHS^P. - ?AMPi c CASE PCP USING 3 ~ ^ 4 

WING 

pcrpp^c c SSOP 

4 

«AL°HA (DEGREES) 



5 

23. 0 



6 

tAFP^rr ^atTO 



7 

1 • 



3 

F T RANFVERS F CU’F 



9 

6. 



11 

C. 1. 2. ?. 

4. 

5. 

1 1 

tSPANH T S- CUTS 



1? 

6 • 



LI 

• 2 .2 • 4 *6 

.9 

1. 

14 

* C r NT FRL T N E 



15 

6. 



16 

*0 r t T A MING PproPOCSSSlR 



17 

C . 



1* 

trprc V0 5 ?T c X 



1 q 

9, 



29 

XDTTCH AXIS 



?1 

C . 25 



?? 

*I*”A T T0N 



?? 

15, 



2 4 

tDpT^ 



r> % 

5* ... j 



?6 

S3 ? c SA C E 
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Beginning with card 16, the input format distinguishes between preprocessors for 
planar delta wing, arrow wing, or gothic wing, and a general cambered wing geometry. 


Card 

Number 

Card 

Column 

Variable 

Name 

(Default Value) 

Description/Comment 

16 

1-25 

or 

or 

$DELTA WING PREPROCESSOR 
$ ARROW WING PREPROCESSOR 
$GOTHIC WING PREPROCESSOR 

17 

1-60 

YLE(I) 

(0.0) 

y-coordinates of wing leading edge 
for each transverse cut; 
this card is intended for use 
with gothic-wing preprocessor 

18 

1-12 


$FREE VORTEX 

19 

1-10 

SFS 

Number of spanwise cuts for 
the free-sheet network 

20 

1-11 


$PITCH AXIS 

21 

1-10 

XPITCH 

x-value of pitch axis 

22 

1-10 


$ITERATION 

23 

1-10 

TMX 

Maximum number of iterations allowed 
for the iteration procedure; the given 
number should be a multiple of 5 
since the AIC’s are updated at 
every 5th iteration 

24 

1-6 


$PRINT 

25 

1-10 

PRINT 

(5.0) 

Printing output occurs at every 
PRINT iteration; the given number 
is recomended to be a multiple of 5 
(see reason given in card 23) 

26 

1-12 


$END OF CASE 


The input format of a general cambered wing geometry is: 


Card 

Number 

Card 

Column 

Variable 

Name 

(Default Value) 

Description/Comment 

16 

1-19 


$INPUT WING NETWORK 

17 

1-10 

FNZ 

Number of corner points; this number 
should be equal to the product of 
number of transverse cuts and 
number of spanwise cuts 
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Variable 


Card 

Card 

Name 


Number 

Column 

(Default Value) 

Description/Comment 

18 

1-60 

ZM(2,I) 

y,z-coordinates of corner points 



ZM(3,I) 

input by each transverse cut 

19 

1-10 

FNLE 

Number of wing leading-edge 
corner points 

20 

1-60 

YLE(I) 

Indices of wing network corner points 
for leading edge from nose to tail 

21 

1-10 

FNTE 

Number of wing trailing-edge 
corner points 

22 

1-60 

YTE(I) 

Indices of wing network corner 

points for trailing edge 

(from centerline to leading edge) 

23 

1-12 


$FREE VORTEX 

24 

1-10 

SFS 

Number of spanwise cuts for 
free-sheet network 

25 

1-11 


$PITCH AXIS 

26 

1-10 

XPITCH 

x-value of pitch axis 

27 

1-10 


$ITERATION 

28 

1-10 

TMX 

Maximum number of iterations allowed 



(5.0) 

for the iterative procedure; 
the given number should be a multiple 
of 5 since AIC’s are updated 
at every 5th iteration 

29 

1-6 


$PRINT 

30 

1-10 

PRINT 

Printing output occurs at every PRINT 



(5.0) 

iteration. The given number is 
recommended to be a multiple of 5 
(see reason given in card 28) 

31 

1-12 


$END OF CASE 


PRINTED OUTPUT DATA 


An example of the printed output format is given in table 3, which shows the results for 
a delta wing of aspect ratio 1 at 20° angle of attack computed in the 15th cycle of the 
iteration procedure. The input data of this case are contained in table 2. The output 
format is self-explanatory, but a few symbols and words must be defined: 

STEP SIZE is the length of the vector Ax (i) obtained from equation (52). 

CIRCULATION ALONG TERMINATED EDGE OF FED SHEET is the strength T of the 
line vortex along the terminated edge computed at the midpoint of the side edge of the 
fed-sheet panel. 
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Table 3 .— Example of Printed Output Data (Continued) 


< *» 

3*£nc .. 

*2121 

.g?3g 

i.:-s<»6 

_ *4_?46 

-.0007 

.0616 

.169* 

- .cro 7 

_,57i5 

-.31*8 

.2527 

. 175C 

10 

’ .5C3C 

.4l?6 

.00?? 

• 99 7 4 

. 9 ; 9<5 

-.0337 

.0479 

.1763 

-• 0 C 37 

1. £72 6 

9226 

. 2 5 . C 

• 1 75C 

1° 

3.6r.jc 

.6105 

. r 7 *5 1 

• °756 

1.1 551 

-.C319 • 

.*171 

. 72 <-2 

-.:r L9 

1.29 6 6 - 

1.365'- 

.2377 

* 175C 

22 

3*5.01113 

... 


l.mi* 

.54*? 

- flfllQ 

^7 7or 

. 646ft 

.P M O 

.61 16 

•.51 67 

* nqqq 

. 1 75C 

21 

4. 5030 

.1125 

• 0030 

1.0064 

.1104 

.0049 

.92C8 

• C 1 1 4 

. 0 [ 1*9 

. 1770 

-.0 250 

.1520 

• 2250 

22 

4.5900 

.3375 

.0 9C0 

• 9961 

• 4 C 84 

• 00 35 

.9123 

• C 4 i 5 

• 00 35 

.3350 

-.1391 

.1659 

.2250 

22 

— h.snnn 

.5425 

.0340 

-9379 

^A?33 

-.P4*9_ 

. O fftl 

-nnn 

• - n r s a 

7(»l63 

• CCjO 

* \ 0 44 

* ??5f 

2% 

4.5030 

• 7075 

•occo 

• 9016 

• 9923 

-.0205 

• 8 49 1 

.3043 

-.0205 

.9940 

-.7979 

• 1861 

• 2250 

25 

4.5000 

1. C 125 

• 0000 

.9754 

.54 33 

-.0108 

.8135 

.511? 

-.0138 

.3234 

-.2468 

.0 767 

• 2250 

26 

Jinan 

-• *1276 

^6656 

i-?06fl 

.41 ?6 

. ???9 

^2?9ft 

.61 65 

. 9 7Q7 

.m ? 


. 


27 

.5000 

• 1306 

.0176 

1.2972 

• 3552 

.3845 

.8128 

.2849 

1.1C14 

.C 321 




2* 

.6000 

.liri 

• 05C1 

1.3150 

• 2586 

.5393 

.0739 

• t 4 C 4 

1.1475 

.0051 




22 

.s flap 

»1?67 

.0421 

l^i412_ 

- 1-1 37 

.6627 

. 914 A 

-.?M 1 

1.1553 

. cn 7? 




1C 

.5000 

• 1 2C T 

.0514 

1.3846 

-• C 768 

.7220 

1.0047 

-.4648 

1.1008 

. 0071 




« 

.5000 

.1122 

.0625 

1.4398 

-.2914 

• 6856 

1.0 857 

-.7375 

. 947 i 

.0350 




- 1? 



-1C16 

-P6ft9 

1.5.:? 5 

• .4925 

.5295 

1.1 755 

- . Qft ^ 1 

.6577 

. P.15Q 




3? 

• 5C3C 

• C 69 ? 

.0714 

1.56 88 

-.6 228 

.2376 

1.2677 

-1.1212 

.1948 

. 0033 




14 

1.6 WO 

.3619 

.317? 

i.22’0 

.3 783 

.2081 

.7505 

.5419 

.9141 

C0C5 




- 1* 

i-5Pao 

- ^2J9S 


i .75:? 

^1243 

. 36P 9 

. 8346 

.2460 

1.P20 7 

-.ring 




IS 

l.rooo 

.3671 

.0911 

1.26*4 

.2 353 

• 5090 

. 0 866 

.-359 

1.P646 

• CJ 26 




*7 

1.5 030 

• 3 766 

.1275 

1.2917 

. 0 99C 

• 6286 

.9 75 9 

-.1024 

1 • C 765 

. CO 56 





l*S0flfl 

- -3586 

.1607 

1^X267 


. 6ft 77 

. 9 94 P 

-.4257 

1 -C 31 3 

.8.1 66 




19 

1.6000 

.1173 

.1682 

1.3744 

-.2031 

.6562 

1.0636 

- • 6 8 C 6 

.8935 

. C369 




40 

1.5030 

. 3tl 3 

.2076 

1.4300 

-.4739 

.5121 

1.1415 

-.5093 

• 63i 2 

• C 3 6 9 




41 

i^5ri£_ 

_ .2*47 

.2156 

1.4925 


.2321 

1 . 2 25 ft 

-1.1479 

.1943 

. ri43 




42 

2.5C30 

.6152 

.029? 

1.1993 

.3 559 

.1887 

.7645 

.4905 

.0813 

-. C019 




41 

2.5030 

• 6464 

• 0900 

1.2092 

.3103 

.3305 

.8313 

. 2314 

.9643 

- . CO 16 




k k 

- 2^5_POtl__ 

.641? 

.1-5? 7 

1.218* 

-A2?l 

.44 31 

-0 746 

. L 3 J ft 

1 .f M 

-.rn n 




45 

2. 5013 

• 622 • 

HRT3 

1.2446 

.0 925 

.5964 

, .9257 

-.1721 

1. P 116 

• 0033 




46 

2.5030 

.5921 


1.2809 

-.0780 

.65 31 

.9842 

-.3978 

.970 3 

. 0051 




... 47 

■VTvTuTS 


wK> H m 


-*£706 

.6244 

PTI731 

-.6351 


B1P4S 4 • I TCTKBI 




44 

2.5030 

.4961 

.1459 

1.3813 

-.4528 

• 4802 

1.1256 

-.8499 


. 0046 




49 

2.5030 

.4351 

• ’50 3 

1.4358 

-.5 74? 

.2236 

1.2C11 

-.98C 4 


. CO 3 2 




so 

__ ’.£*30 

.5929 

-1391 

1.1452 

.3731 

.1644 

,.7677 

.5497 

HViJ J|l 1 pi 

■mim 




51 

1.5030 

fnps^ 

.1226 

1.1059 

.3289 

.3474 

.7774 

.2741 

.6623 

C061 




5? 

’.5800 

BBC tQjU 

• 2106 

1.1124 

. 2 5 OC 

• 4860 

.5251 

• C666 

.9498 

-.0394 




53 


HR 


1-1542 _ 

.1216 

.6112 

-47fl 

-.1174 

.9777 





54 

3.5030 

. *449 

• ’77 1 

1.2336 

HTI 

.6599 

.9731 

mi 

.9497 

-.0048 




55 

’.5000 

.7<61 

.417* 

1.2652 


.6332 

1.012? 

mBt 3 

.0345 

to : 6 




__ 56 

wmvwwm 

..*1112. 

■HlXtl 

■VnuS 

WORK t SI 

.4974 

■Klim 

■LIVOi K 


IWJiTXKI 




57 

3.5C30 

.6264 

.50’3 

1.3859 

-.5425 

.2361 

1.1711 

-.9255 

. ?19 3 

. c: 3; 




5* 

4.6rao 

1.140“ 

.0499 

1.-618 

.3 766 

.20 61 

.5266 

.516? 

.*C?4 

- . CO u 9 




50 

4.5030 _ 

_ 1.162 7 _ 

.1559 

1.2904 

.’409 

.3691 

. 8 ’24 

__*£4 Ofl 

.7492 





SC 

4.6C30 

1.1859 

■KTTl 

mm 

WKBM 

• 5201 

.*491 

.1704 

mmum 

.0032 




61 

4 . r C3 C 

1.16»7 




.6290 

.5 59i 

-.:236 


- • C3 1 c 




62 

. k -5_c: :■ 

1 .1 1 f 2 



^wn m 


.* cr 7 


■mm 

j £,5 




61 

4 . c 03 P 

t . 0 4 4 ? 

.568’ 

1.1492 

-. 2384 

• 6463 

.9475 

- • 5 C 0 2 

.8077 

- • CO 75 




64 

4.6000 

• ?5f 3 

.5**03 

1.2131 

-.4167 

.5095 

1.0198 

-.7177 

.5940 

-. t041 




J65 

4.5C0Q 

i £A 1 2 - 

.-^574 

U2A IB 

-*1227 

.25 36 

-1*1022 

,-.^3S63L 

-• 2.422— 






TC T 5 NT S • T 965 

cj r r hing mgmfnt co'T^inirNt = .42?e 

pirx x_rs _5 : : *25 ax 

COOT CHC-0 = 5.occo 

WTHG = 6 • 25 C 0 









Table 3.— Example of Printed Output Data (Concluded) 


X Y Z COORHINATc re CCRNEQ POINTS 


FRFE SH?" N^ThORK 


. 000 

.700 

. ooc 

. 700 

• 90 7 

• 007 

.000 

.000 

• 00 0 

.COC 

.ooc 

.000 

• OOC 

.00 u 

• OCC 

. ooc 

.000 

.-300 

. 70 P 

• 79 C 

• 00 0 

.CC-0 

.000 

.00 0 

• COO 

.occ 

.000 

l.oor 

.250 

.0 00 

IaQDC 

*-2611 

^323._ 

■ --1.2 Cl 

.262 

_ ..0 4ft 

1-f c ' 

.250 

.1/72 

1 . t_ lj ’ 

. 24© 

. " 96 

! « ' L 

. 2 3u 

.117 

1.9C0 

.215 

. 133 

1. OCC 

.191 

.14? 

1.C0C 

.166 

. 143 

2.00 0 

• 5 0 C 

.000 

2.0CZ 

.517 

.046 

?. OOP 

. 51 9 

. 397 

2. 700 

.59® 

.147 

2. CCO 

.490 

. 194 

2 . C C C 

.461 

.235 

2. OOC 

.423 

.268 

?. non 

* ??6 

.-232 

? . j r r. 

. *2 5 

.297 

3.LC P 

.75 0 

.rn r 

3. COC 

.774 

. 0 7 n 


' § Z ?6 

. 1 uf> 

3. 700 

. 761 

. 221 

3. PCC- 

.“>31 

.291 

3.000 

.696 

. 352 

3.wC: 

.620 

. 40 : 



3 . c 

.5 53 

.429 

3. OOC 

. 4*1 

.437 

u. or: 

l.n^f 

' .PC 7 

4.CC 7 

1.749 

.707 

4.PCC 

1.050 

.107 

4. COC 

1.045 

• 280 

u _ or * 

-1^10 

. 35 3 

,4 a coc . 



.-467 

4.00 0 

.0 77 

. 532 

4 . 0 0 Q 

. 704 

.5 73 

4 . 0 c r 

^£32 

. 50C 

5. 0PC 

1. 250 

• OOP 

5. 7GC 

1.3C1 

.113 

5.CCC 

1.324 

.237 

5 . C- 0 0 

1.317 

.364 

5. COD 

1.283 

.406 

5. 00G 

l v . ?1« 

.595 

5.7C7 

1.127 

.651 

5.CC C 

1.012 

. 73 4 

5.0Cc 

. 005 

.742 











Ff C SHE£ T 

NE 7WC p K 








JXC 

A jn r 

.cc: 

. M r 


1.CC3 

.166 

.14 3 

1. GOG 

.162 

.G 79 

2 . aor 

.325 

.207 

2.700 

.m 

.159 

2. OOC 

.491 

.427 

3. COO 

.471 

.238 

*.coc 

.68? 

• 5 SC 

4 . 00 c 

.669 

.324 

5. OOC 

• 9*5 

. T 42 

5.:tc 

.56* 

.422 










CIRCULATION ALONG WING TRAILING EDGE is computed at the midpoints of the 
panel trailing edges. 

zcx 

= x 0 ) 



ZCY 

* yo 

► 

coordinates of panel center 

zcz 

= z 0 ) 



vux 

- (V »>x 1 


components of V u on upper side of the 

VUY 

= CV U )y 

► 

sheet surface in x,y ^-coordinates 

vuz 

= (V«), 1 



VLX 

- (V l } x 1 


components of Vj on lower side of the 

VLY 

= (Vi) 

► 

sheet surface in x,y,z-coordinates 

VLZ 

- (V 1 } Z J 



DCP 

= Acp 


jump in pressure coefficient across the sheet 

CPU 

= C P U 


pressure coefficient on upper wing surface 

CPL 

= Cp l 
= Sj 


pressure coefficient on lower wing surface 

AREA 


wing panel area 
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VERIFICATION OF THE METHOD 


RESULTS 

Numerous example cases have been executed to validate the method and its generality. 
Cases selected are compared with available theoretical and experimental data for a 
range of different geometric configurations including delta, gothic, and arrow wings. 

The capability of the method to predict overall wing coefficients is shown in figure 20 
for a delta wing of aspect ratio 1 at low subsonic speed (Moo = 0). This figure shows the 
well-known nonlinear variation of the normal-force coefficient cn with angle of attack. 
Several values of cn were computed for angles of attack up to 20° and agree well with 
the experimental data of Peckham (ref. 17) and theoretical results from the 
leading-edge-suction analogy of Polhamus (ref. 3). The corresponding load distribution 
at a = 20° is plotted in figure 21 and compared with Peckham’s experimental results. 
Although only 25 wing panels were used on one-half of the configuration, the 
completely three-dimensional, nonconical load distribution was well predicted, including 
the location of the vortex-induced pressure peaks and the decrease of the load toward 
the trailing edge. 

Figures 22, 23, and 24 show detailed surface pressure distributions for a delta wing of 
aspect ratio 1.4559 at a = 8.8°, a = 14°, and a = 19.1°. Upper and lower surface pressures 
are well predicted for the higher angles of attack, as the comparison with experimental 
data (ref. 18) illustrates. At 8.8° the differences, although unclear, may be due to the 
blunt trailing edge of the experimental model or an inadequate definition for the 
fed-sheet geometry. The experimental results clearly show the effect of the secondary 
vortex separation, which takes place on the upper surface just slightly outboard of the 
main vortex. The present method does not model secondary vortex separation and, 
consequently, produces a slightly different shape for the pressure peaks. 

The method can be applied to more general configurations. For example, figure 25 
shows the method applied to a gothic wing having a swept trailing edge and a curved 
leading edge. This figure shows good agreement of the normal-force coefficient cn with 
experiments (ref. 17) at the relatively high angle of attack of 14.3°. 

Figure 26 shows the predicted load distribution of an arrow wing. Experimental data 
are not available for comparison, but the plotted loads appear to be realistic and 
demonstrate that the method is capable of handling other than delta-wing planforms. 

CONVERGENCE CHARACTERISTICS 

The progress of the solution is monitored by examining the residual errors in the zero 
pressure jump boundary condition on the free sheet and wake, and the stream-surface 
boundary condition on the wing. Figure 27 shows the convergence characteristics for a 
delta wing of aspect ratio 1 and a gothic wing of aspect ratio 1.60. Each configuration 
has 25 wing panels. 
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Of, deg 

Figure 20. —Normal Force Coefficient C/g, Delta Wing 
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AR = 1.4559, M co = 0, a = 14.0 


C P k 


— — Present method 
O A O □ Experiment of 
Marsden 



Figure 23.— Surface Pressure Distribution of Delta Wing, a =14.0° 
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AR = 1.4559, = 0, a = 19.1° 

" Present method 

O A O O Experiment of 
Marsden 

30 Wing panels 



Figure 24.— Surface Pressure Distribution of Delta Wing, a = 19. 1° 
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a = 15.8° M m = 0 

= 71.2° c r = 111.96 cm 

b = 112.88 cm 






•25 h 


AR = 1 .60, = 0, A te = 56.3° 

25 Wing panels 


0 


± 

5 


_L 

10 




Number of iterations 



.25 


a = 14.3° 

AR = 1, = 0 

25 Wing panels 


I i ■ I I I I fr, 

0 1 2 3 4 5 6 

Number of iterations 

Figure 27. -Convergence Characteristics 
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The paneling and convergence characteristics for the delta wing are shown at the lower 
side of figure 27, where the normal-force coefficient cn is shown as a function of 
iteration number for two different angles of attack. The dashed lines indicate the value 
of cn obtained from the leading-edge-suction analogy. For this case, the aerodynamic 
influence coefficients were updated in each cycle of the iteration. The solution quickly 
seeks a level after only one or two iterations and then exhibits some oscillation, often a 
characteristic of Newton schemes. 

The paneling and convergence characteristics for the gothic wing are shown at the 
upper side of figure 27. For this case, the aerodynamic influence coefficients (whose 
computation consumes the largest fraction of computer run time) were updated only 
after the fifth and tenth iterations. The difference in convergence characteristics for the 
two cases is apparent. 
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CONCLUSIONS 


The work reported here demonstrates one of the applications of a new general, subsonic 
potential flow computational technique recently developed (reported in ref. 14). With 
the use of this technique, a three-dimensional method has been formulated for 
predicting the flow field about swept, sharp-edged wings characterized by the presence 
of vortex separation at the leading edge. The basic approach has been verified for 
selected wing planforms such as delta wings of different aspect ratios, a gothic wing, 
and an arrow wing. Further numerical experiments conducted after completion of the 
contract work uncovered convergence problems with the present computer code for more 
general wing geometries like double delta-wing planform, diamond wings, and cropped 
delta wings. The problems do not appear to be fundamental in nature and are due in 
part to a lack of experience concerning panel density requirements, the accuracy of the 
initial guess, and the appearance of a second vortex at a planform break that is not 
modeled in the present code. The work to date has been successful in overcoming the 
most difficult aspects of the problem and provides the basis for development of this 
initial capability into a method suitable for supersonic flow and for complete 
configuration analysis in subsonic and supersonic flow. In addition, it can be applied to 
related problems involving free vortex sheets that are encountered in the areas of 
powered lift, jet interactions, jet flaps, and so forth. 


Boeing Commercial Airplane Company 
P.O. Box 3707 

Seattle, Washington 98124 
September 24, 1975 



APPENDIX A 

COMMENTS ON CHOICE OF FED-SHEET MODEL 


A simple method for treating the core of the rolled-up primary vortex sheet was used by 
Smith (ref. 6) in his conical approach to the leading-edge vortex problem. This model 
consists of a single vortex with a fed sheet (see fig. 28), which he calls "feeder sheet.” 
The sheet size and orientation are adjusted so that the total force on the vortex and on 
the fed sheet is zero to provide a unique solution. 

The three-dimensional method described in this document uses an even simpler 
representation of the viscous core region. The fed sheet is treated as an entirely 
kinematic extension of the free sheet, and no boundary condition is applied to the fed 
sheet. The tacit assumption in this simplified model is that boundary conditions applied 
to the free sheet are sufficient to adequately position the fed sheet, whose typical 
dimension is small compared with dimensions associated with the free sheet and with 
-distance from the wing. Size and initial position of the fed sheet are taken from conical 
flow results of Smith. 

This fed-sheet model was the outcome of a theoretical investigation of the relation of 
mathematical discrete line vortices to the problem of leading-edge vortex separation. In 
particular, it was discovered that Smith’s fed-sheet model, when applied to the 
three-dimensional nonconical flow problem, would contain two self-induced infinite 
forces, with no possibility of mutual cancellation, and that for this reason the boundary 
condition of zero integrated force over the vorticity could not be applied (ref. 19). 

The origin of these infinite forces can be explained with the aid of figure 29, in which 
the fed sheet carrying a surface distribution of vorticity y is shown. This vorticity 
merges into a discrete line vortex of strength T along the terminated edge. The line 
vortex along that edge will be curved in the three-dimensional and nonconical flow 
problem. This curvature is responsible for one of the infinite forces, since any curved 
discrete line vortex induces on itself an infinite component of velocity. When a curved 
vortex is free to move with the surrounding fluid (as is the case with free leading-edge 
vortices), it will move with infinite speed and will change, in general, its shape with 
infinite speed (ref. 19). The curved discrete vortex forming the terminated edge is not 
allowed to move with infinite speed and, for that reason, will always sustain a trans- 
verse component of force that is infinite in magnitude. 

The other infinite force is experienced by the fed sheet and acts in the direction parallel 
to the terminated edge. It is present in all cases where the fed sheet is feeding vorticity 
into the edge vortex, i.e., when the strength of the edge vortex changes along its length. 
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APPENDIX B 

PANEL INFLUENCE COEFFICIENTS 


B.l INTRODUCTION 

In this appendix we shall calculate the potential and velocity induced by a source or 
doublet distribution on a curved panel. As shown in figure 30, let S be the curved panel 
surface, 2 its tangent plane projection, Q a point on S, n the normal to S at Q, and ? a 
field point. The potential <f> at ? induced by a source distribution a on S is defined by 

♦■JOT 

s 

where 

R = (£ - *,T7 

and 

R = IRI = V(i - x) 2 + (t? r y) 2 + (f - z) 2 

The potential (/> at P induced by a doublet distribution /x on S is defined by 

-4 

. =//s 

The velocity induced at P by a source or doublet distribution on S is defined by 


MVq 


4ttR 


•n dS 




i? A ' 

-R • n 
4 ttR 2 


dS 


(B-3) 


4irR 


dS 


w y m d 


(B-l) 


m o\ 


v = V (j> 


(B-4) 


We assume that the surface S is defined by 

f = a + b 17 “ , . (B-5) 

We also assume that S does not deviate significantly from 2, more precisely that 

5 << 1 (B- 6 ) 


where 

r l Max 5 \\ T3 ~7j? Jl. , n 

8 = 2(£,,)eZ S- + b-,-( . 

(Nominally we assume S < 0.066.) 
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The distribution of singularity strength on S is assumed to be linear in the case of a 
source panel and quadratic in the case of a doublet panel. To be specific we assume 

a = o Q + £ + a v t] , Q, V ) e 2 (B-8) 

and 

l 2 1 j 

M = Mo + MgS + + + J M Vr} V > (€,*?) e S . (B-9) 

B.2 EVALUATION OF SOURCE AND DOUBLET 
INTEGRALS FOR AN ARBITRARY FIELD POINT 

Let us first consider the evaluation of source potential defined py equation (B-l). 
Evaluation of source velocity and doublet potential and velocity will quit^e similar. 
The first step in the evaluation procedure is to transfer the integral over S to the 

(BrlO) 


(B-ll) 

Substitution of equations (B-2), (B-8), and (B-ll) into (B-10) yields an explicit integral 
for cf > . However, the integral cannot be evaluated in closed form as it stands. By 
employing the hypothesis that 8 2 is negligible compared to unit (hypothesis (^-6)) the 
integrand can be approximated by terms that are integrable in closed form. A qniform 
approximation to sec (£,n) can be obtained by noting that 

4a 2 £ 2 + 4b 2 r? 2 <1651 (B-12) 


equivalent integral over 2. We have 

/ -1 \ AA 

</> = ff °\ 4 ^) sec d $ dT ? f 
2 

From equation (B-5) we obtain 

sec (f,n) = Vl + + 4b^??“ . 


hence 


AA 

sec (f ,n) 


(B-13) 


A uniform approximation to 1/R is somewhat more difficult to obtain sin^e this factor is 
singular. Let (xo,yo) be the point on 2 closest to (x,y) and set 


z o = ^ x o> y o^ h = z-z o and r = # - x ) 2 + (v - yfi 


(B-14) 


Then 


R = V(r 2 + h 2 ) - 2h(f - z 0 ) + (f - z 0 )^ 



1 + 


-2h(f - Zq) + (i‘ r Zq)- 
r 2 + 


(B-15) 
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Let 


e = 


0‘ 


Then 


M a x \ j£ z i 

(£, r?) e E | r J ' 

-2h(f - zq) + (? - zq)^ I 2erh + ^ ? 

^ r ^ € + € 


r 2 + 1,2 


r 2 + 1, 2 


Therefore if e 2 is everywhere negligible, equation (B-15) yields 

11. - z o> •• 


R p + 


= . 


But 


Max ’ jlf(^)-r(x 0 ,y 0 ) 


Max j I a (£ + x 0^ - x o) + + yo^ 7 ? ~ yo') l 


(B-16) 


(B-17) 


(B-18) 




sC 


V a2 (^ + x q)" + b 2 (r? + yo) 2 


Max 

(£>f?) e 2 

\ 

|Va 2 ({ + x 0 )2 tb 2(, + yo) 2 . 2 j 


w - x o )2 + 0? - yo>" 


(B-19) 


< 85 . 


A much better bound on e is available when (x,y) is several panel diameters away from 2 and 
the assumption that 8 2 is negligible becomes unnecessary. However, in this case a far-field 
expansion will be used to obtain an efficient approximation to the right side of equation 
(B-10). 

Substituting equations (B-8), (B-13), and (B-18) into (B-10) and rearranging we obtain 

4> = o(x, y) 1(1, l) + o x (x, y) 1(2, 1) + a y (x, y) 1(1,2) (B-20) 


where 


o(x, y) = Oq + a^x + a^y (B-21) 

o x (x, y) = a£ 

ffy( x ,y) = o n 
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Here 


I(M,N) = -^t[h(M,N, 1) + a(hH(M + 2, N, 3) + 2xhH(M + 1 , N, 3)^ 

+ b(hH(M, N + 2, 3) + 2yhH(M, N + J , 3)) (B-22) 

+ c(hH(M, N, 3))J . 

where 

c = a x^ + b y~ t zq (B-23) 


and 


H(M, N, K) 



(£-x) M 1 (r? - y) N 1 I} . ^ 

/ F= - 5 tW d S dr ? 

(V(l - x) 2 + (t? - y) 2 + h~J 


(B-24) 


The H integrals will be evaluated in the next section. 

To find V, equation (B-20) can be differentiated. For this purpose zq may be treated as a 
constant since its derivatives with respect to x and y either cancel or are negligible on 
account of hypothesis (B-6). The derivatives of the H integrals, fhep, are simple 
combinations of the H integrals themselves; i.e., 

N, K) - -(M- 1) H(M - 1, N, K) + KH(M + l.N, K + 2) 


H(M, N, K) = -(N- 1) H(M, N- 1, K) + KH(M, N+.l, K + 2) (B-2$) 


~ H(M, N, K) = -KhH(M, N, K + 2) . 

Actually, it turns out to be easier to calculate V by differentiating equation (B-10) and 
using a generalized form of equation (B-18); that is, 


j \_ Kli(C-z 0 ) ' 

R K %K + „K« • " ' 


- P 


^ + h- . 


In either case we obtain 

V = o(x, y) T( 1 , l) + o x (x,y)?(2, 1) * o y (x. y)T(l, 2) 

where 

7(M, N) = (j x (M,N), J y (M,N), J Z (M,N)) 


(B-£6) 


(B,27) 


(B-?8) 
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and 


; J X (M,N) = -^[H(M + 1,N,3) + a(3hH(M + 3, N, 5) + 6xhH(M + 2 J .N,;5)) 

+ b(3hH(M+ l,N,+ 2,5) + 6yhH(M+ 1,N+ 1,5)) 
+ c(3hH(M+ 1,.N, 5))] 

J y (M,N) = -^[h(M,N+ 1,3) + a(3hH(M + 2,N+ 1, 5) + 6xhH(M + 1,N + 1, 5)) 

+ b(3hH(M, N + 3, 5) + 6yliH(M, N + 2, 5)) 

+ c(3hH(M, N+ 1, 5))J 

J Z (M, N) = - 4^ [~ hH(M, N, 3) + a(H(M + 2, N, 3) - 3h 2 H(M + 2, N, 5) 


+ 2xH(M + 1 , N, 3) - 6xh 2 H(M + 1 , N, 5)) 
+ b(H(M, N + 2, 3) - 3h 2 H(M, N + 2, 5) 
+2yH(M, N + 1 , 3) - 6yh 2 H(M, N + 1 , 5)) 
+ c(h(M, N, 3) - 3h 2 H(M, N, 5))] . 


Doublet potential and velocity can be evaluated similarly starting with equations (B-3) 
and using 

n = r 1 (- 2a£, - 2br?, l) (B-29) 

yi + 4a 2 £ 2 + 4b 2 7? 2 

along with equation (B-26) and hypothesis (B-6). We obtain 


<(> = 'ju(x, y) 1(1, l) + ju x (x, y) 1(2, 1) + M y (x,y) i(l,2 ) 

+-^H xx (x,y) 1(3, 1) + M xy (x, y) 1(2, 2) +-TM yy (x, y) 1(1,3) 

where 

1 0 1 9 

p(x, y) = M 0 + ju^x + n v y + -M^x- + p^xy ^-n^y- 
M x (x, y) = Hj; + M||X + M|^y 

My (x,y) = ’ • 

. M xx (x, y) = n# 

M xy (?<, y) - 

Myy(X,y) ; = 


(B-30) 


(B-31) 
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and 

KM, N) = 4^r[hH(M, N, 3) i a(H(M + 2, N, 3) + 3h 2 H(M + 2, N, 5) + 6xh 2 H(M + 1, N, 5)) 

+ b(H(M, N + 2, 3) + 3h 2 H(M, N + 2, 5) + 6yh?H(M, N + 1 , 5)) 

’ (B-32) 

+ c(-H(M, N, 3) + 3h 2 H(M, N, 5))] . 

We also obtain 

V = M(x,y)T(l, 1) + M x (x,y) 1(2, 1) + M y (x,y)7(l,2) (B-33) 

+3M xx (x, V) 7(3, 1 ) + M xy (x, y) 7(2, 2) +4^yyf x - y)?(l, 3) 

where 

7(M, N) = (j x (M, N), J y (M, N), J ? (M, N >) (Br34) 

and 

J X (M,N) = ^jr[3hH(M+ 1, N, 5) + a(3H(M + 3, N, 5) - 2H(M + l, N, 3) 

+ 1 5h 2 H(M + 3, N, 7) - 2xH(M, N, 3) + 30 xh 2 H(M + 2, N, 7)) 

+ b(3H(M + 1 , N + 2, 5) + 1 5h 2 H(M + 1 , N + 2, 7) + 30yh 2 H(M + 1 r N f f , 7)) 
+ c(-3H(M + 1, N, 5) + 15h 2 H(M+ 1,N,7))1 

J y (M, N) = ^-[3hH(M,N + 1,5) + a(3H(M + 2, N + 1,5)+ 15h 2 H(M + 2,N + 1,7) 

+ 30xh 2 H(M + 1, N + l,7)) + b(3H(M.N + 3,5)-2H(M,N + 1,3) 

+ 1 5h 2 H(M, N + 3, 7) - 2yH(M, N, 3) + 30yh 2 H(M, N + 2, 7)) 

+ c(-3H(M, N + 1 , 5) + 1 51 i 2 H(M, N + 1 , 7))] 

J ? (M, N) = ^[h(M, N, 3) - 3h 2 H(M, N, 5) + a(3hH(M + 2, N. 5) 

- 1 5h 3 H(M + 2, N, 7) + 1 2xhH(M + ) , N, 3) - 30xh 3 H(M + 1 , N, 7)) 

+ b(3hH(M, N + 2, 5)- 1 5h 3 H(M, N + 2, 7) + 12yhH(M,N+ 1,5) 
t 30yh 3 H(M , N + 1,7)) +c(9hH(M, N, 5)- 15h 3 H(M, N, 7))] . 
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B.3 CALCULATION OF H INTEGRALS 


In this section we shall compute in closed form the integrals 

H(M, N, K) = ff ^~ X)M ~ 1 ( J?-y> N-1 d| dr?, 

p = - x ) 2 + (t? - y) 2 + h 2 

for M = 1, MXQ; N = 1, MXQ-M + 1; K = 1, MXK,2. 


(B-35) 


Here MXQ is the maximum value of M + N-l required and MXK is the maximum value 
of K required. The following values of MXQ and MXK are evident from the previous 
section. 


Panel Type 


MXQ MXK 


Source potential (flat panel) 
Source velocity (flat panel) 
Doublet potential (flat panel) 
Doublet velocity (flat panel) 
Source potential (curved panel) 
Source velocity (curved panel) 
Doublet potential (curved pianel) 
Doublet velocity (curved panel) 


2 

3 

3 

4 

4 

5 

5 

6 


1 

3 

3 

5 

3 

5 

5 

7 


(B-36) 


The integrals H(M,N,K) may be computed with the aid of the following algebraic 
recursion relations. We have the obvious identity 

H(M + 2, N, K) + H(M, N + 2, K) + h 2 H(M, N, K) = H(M, N, K - 2) (B-37) 


Integration by parts yields 


(K - 2) H(M, N, K) = (M - 2) H(M - 2, N, K) - £ ^ F(M - 1 , N. K - 2) (B-38) 


and 


(K - 2) H(M, N, K) = (N - 2) H(M, N -2, K - 2) - £ v F(M, N - 1 , K - 2) . (B-39) 


The summations on the right sides of equations (B-38) and (B-39) are over all four sides 
of 2 with the contribution of a typical side L displayed. Here v is the unit outer normal 
of the side L (see fig. 31) and F(M,N,K) is the line integral defined by 


F(M, N, K) - 


c (| - x) M_I (7? - y) N_1 

AO , n - \l 

' n K 

ux 1 fj — y 


- V(S-x ) 2 + (r?-y ) 2 + h 2 • (B-40) 
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Figure 31 .-Quadrilateral Geometry 
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The procedure for evaluating the F integrals will be described later. Assuming the F 
integrals are known, the recursion relations (B-37), (B-38), and (B-39) may be 
recombined to yield an efficient procedure for calculating the H integrals. Because some 
of the H integrals are singular on the edge of S, it is actually necessary to consider 
three slightly different procedures, depending on the relationship of the Field point to 
the panel. Let us define d^ to be the minimum distance of the point (x,y) to the 
perimeter of 2. If 8h is some small number (nominally chosen as 0.01), we have the 
following three procedures. 

Procedure 1: |li[ > 5j a djj 

1. 4 4 

H(l, 1, 1) = — l h I £ tan-^a^C! -fi 1 c 2 ) ! c 1 c 2 + a 2 £ 1 £ 2 l + £ aF( 1 , 1,1) (B-41) 
1 1 

where 

cj = g 2 + |h| V^ 2 + g 2 , c 2 = g 2 + |h| yC^ + g 2 

and 

g = V a 2 + h 2 


2. 


H(l. 1, K) = 


1 


(K - 2) h“ 


(K - 4) H( 1 , 1 , K - 2) + £ aF(l, 1 , K - 2) 

1 


K = 3, MXK, 2 . 


(B-42) 


3. 


H(2,N, 1) = 


1 


(N + 1) 

N = 1 , MXQ - 1 . 


h 2 £ ^F(1,N, 1)+ £ aF(2, N, 1) 
1 1 


(B-43) 


4. 


H(l, N, 1) = 


-h 2 (N - 2) H( 1 , N - 2, 1) + h 2 £ i^F(l,N-l,l) 

1 


+ £ aF(l , N, 1) 
1 


; N = 2, MXQ . 


(B-44) 
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5. 


H (M, N, 1) = ( ' m + N - 1) 


-h 2 (M - 2) H(M - 2, N, l) + h 2 £ »^F(M - 1, N, 1) 

1 


+ £ aF(M 5 N, 1) 
1 


(B-45) 


; M = 3, MXQ and N = 1.MXQ-M+ 1 . 


6 . 


H(1 , N, K) = 


1 


(K.-2) 

N = 2, MXQ and K = 3, MXK, 2. 


(N - 2) H( 1 , N - 2, K -2) - £ ^F(l, N- 1, K-2) 


(B-46) 


7. 


H(2,N,K) 


1 

(K-2) 



^F(l,N,K-2) ; 


N = 1, MXQ -1 and K = 3, MXK, 2. 

H(M, N, K) = -H(M - 2, N + 2, K) - h 2 H(M - 2, N, K) + H(M - 2, N, K - 2); 
M = 3, MXQ and N = 1 , MXQ - M + 1 and K = 3, MXK, 2. 


Procedure 2: |h| < S^djj and (x, y) e 2 

1. H(l, 1 , NHK + MXK) = 0.0 

where NHK is a positive integer (nominally taken to be 16). 

2 . 

H( 1 , l.K-2) = YkT4) 

K = NHK + MXK, 3, 2. 

3. All remaining steps are as for Procedure 1. 

Procedure 3: |h| < S^dp^ and (x, y) e 2 
Define 

H*(M, N, K) = H(M, N, K) - 2ttX(M, N, K) |h| M+N " K 


h 2 (K - 2) H(1 , 1, K) - £ aF 0> l.K-2) 

1 


(B-47) 


(B-48) 


(B-49) 


(B-50) 


(B-51) 
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where 


X(M, N, K) 


0 if M or N is even 


[l • 1 • 3 - 5 . . . (M - 2)] [l • 1 • 3 • 5 . . . (N - 2)1 
(K - 2)(K -4)(K - 6) . . . (K - M - N)] 


otherwise. 


Then Procedure 2 may be used to calculate H* . 

We now evaluate the integrals F(M,N,K) for the indices M,N,K required by the H 
evaluation procedures. It is apparent that we need only the following F integrals: 

F( l, 1, K) ; K = 1, MXFK, 2 ' ' 


F(M, N, f) ; M = 1 , MXQ and N = 1 , MXQ - M + I 


(B-53) 


where 


F(1,N,K) ; N = 2, MXQ and K = 3, MXK - 2, 2 


MXFK = 


MXK -2 if |h I > 5 h d H 
NHK + MXK-2 if |h| < § h d H 


(B-54) 


These integrals may be obtained with the aid of three recursion relations. We have the 
following obvious identities: 

F(M + 2, N, K) + F(M, N + 2, K) + h 2 F(M, N, K) = F(M, N, K - 2) (B-55) 


^F(M + 1 , N, K) + ^F(M, N + 1 , K) = aF(M, N, K) . 


(B-56) 


Integration by parts yields 


-(M - l)i»F(M - 1, N, K - 2) + (N - 1) ^F(M, N - 1, K - 2) 

HK-2)u F(M+ 1,N,K)-(K-2)^F(M,N+ l.K) = E(M,N,K-2) 


(B-57) 


where 


E(M, N, K) = 


(£ - x) M 1 (r? - y) N 1 
«K 


(B-58) 


p = V (? - x ) 2 + (*? - y ) 2 + h 2 • 
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The quantities E(M,N,K) may be evaluated directly or else recursively with the aid of 
the formula 


P(I) = (x 2 + Xj)P(I- 1 )-x ] x 2 P(I-2) 


(B-59) 


where 


P(I) = A 2 x 2 * 1 - AjXj 1 


The recursion relations (B-55), (B-56), and (B-57) may be recpmbined to yield an 
efficient procedure for evaluating the required F integrals. Here, again, the singular 
behavior of some of the F integrals near the edges of S requires a special case. Let us 
define dp to be the minimum distance of the point (x,y,z) to the perimeter of S. TJien if 
8 g is some small number (nominally chosen as 0.01), we have the following two 
procedures. 


Procedure 4: g > 5gdp 

1 . 


F(l, 1, 1) = logl v* z +g z + e, 


(v4 2 + g 2 + #) 


= < 


/y £ 2 2 + g 2+£ 9 

logl - /» «j.«2 > 0 


/V«i 2 + g 2 - 8 i 

log ]==f -|; «1,*2 < 0 

W + g 2 - fi 2 / 

(v c i 2 + g 2 - 6 l)(V fi ? 2 + g 2 +«2 


log' 


4 


(Br60) 


i 2 > Q.i\< 0 


2 . 


1 


F( 1 , 1 , K) = — 

g 2 (K -2) 

K = 3, MXFK, 2. 


[(K-3)F(1, l,K-2)-^E(2, 1, K - 2) + ^E(l , 2 , ^ - 2)] ; 


(B-61) 
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(B-62) 


3. a. If 1^1 < |^| 

(0 F(l, N, 1) = (jq-r T y[(2N-3)a^F(l 5 N-l, 1) 

,-(N-2)(a 2 + ^ 2 h 2 ) F(l,N-2, 1) + ^E(1,N- 1, 


N = 2, MXQ. 

(ii) F(M, N, I) = F(M - 1, N + 1, 1) +^F(M - 1, N, 1); 

v t V $ 

M = 2, MXQ and N = 1 , MXQ - M + 1 . 
b. If 1^1 < 1^1 

(i) F(M, 1,1) = (mTT) [(2M -3) a^F(M -1,1,1) 

-(M-2)(a 2 + ^ 2 h 2 )F(M-2, 1, 1)-^E(M- 1, 1 

M = 2, MXQ. 


(ii) F(M, N, 1) F(M + 1 , N - 1 , 1) +—■ F(M, N - 1 , 1); 
T? V 

N = 2, MXQ and M = 1 , MXQ - N + h 


F(1,2,K) = iy»F(l,l,K) « K ~ T > 

(K-2)(6 2 +g 2 )—f~ , 

K = 3, MXK - 2, 2. 

5. 

F(1,N, K) = 2a^F(l,N- l,K)-(a 2 + ^ 2 h 2 )F(l,N-2, K) 

+ ^ 2 F(l,N-2,K-2) ; 

N = 3, MXQ and K = 3, MXK -2, 2. 

Procedure 5: g < 5gdp 
1. F(l, 1, MXFK + NFK) = 0 


where NFK is a positive integer (nominally taken as 16) . 


»]•, 


(B-63) 


1 )]; 

(B-64) 


(B-65) 


(B-66) 


(B-67) 


(B-68) 
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F( 1 , 1 , K - 2) = (k 1 3) [g 2 (K-2)F(l, 1,K) + ^E(2 S 1, K - 2) - ^E(l, 2, K - 2)} ; 

K= MXFK + NFK, 5, 2. (B ' 69) 

3. F(l,l,l) as well as all other F integrals may be computed in the same manner as 
for Procedure 4. 

This completes the calculation of the H integrals. 


B.4 EVALUATION OF SOURCE AND DOUBLET 
INTEGRALS FOR A DISTANT FIELD POINT 

If the point P is a large distance from S, the approximation (B-18) may be replaced by 
an approximation based on this fact. Let 

P = I Pi and Q = |(^| (B-70) 


Then 


Let 


Then 


Hence, if 



e << 


(B-71) 


(B-72) 


(B-73) 


(B-74) 


we have 

r_ i .(?•$) i/q 2 3#* q) 2 

R ~ P p 3 2^ p 3 p 5 

Only the first three terms of the binomial expansion are displayed (monopole, dipole, 
and quadrupole). In practice this expansion is used only when e is less than 1/5. All 
three terms are used unless € is less than 1/8, in which case only the first two terms are 
required. 

Substituting (B-75) into (B-IO) and using hypothesis (B-6), we obtain for the source 
potential 

0 = a 0 I(l, 1) + ajl(2, 1) + o v l{l, 2) (B-76) 
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Here 


where 


and 


I(M, N) = 


C(M, N) + (q(M, N) v^) 

P 2 


1 1 


-~(qq(M,N)-3p • pq(M 


,N)) 


C(M, N) = JJ t?*" 1 d| dr? , 

‘IC 


q(M, N) = (C(M + 1, N), C(M, N + 1), aC(M + 2, N) + bC(M, N + 2)) , 

qq(M, N) = C(M + 2, N) + C(M, N + 2), 

pq(M, N) = (p - q(M + 1, N), p'*'q(M, N + 1), ap • q(M + 2, N) 

. + bp • q(M, N + 2)), . 


_ P 
p P 


(B-77) 


(B-78) 

(B-79) 

B-80) 

(B-81) 

(B-82) 


In computing the last component of pq, the last component of may be ignored on 
account of hypothesis (B-6). The integrals C(M,N) will be evaluated later. The induced 
source velocity may be obtained by differentiating equation (B-76) and is given by 


V = a 0 7(l, 1) + \) + o v T(\, 2) 

where 

. J(M, N) = ^ C(M, N) P + ^(q(M, N) - 3(?^(M,.Np) 

+ -p r (3p^(M, N) +|(qq(M, N) - sfc • p^(M, N)))?) . 


(B-83) 


(B-84) 


A similar expansion may be obtained for the doublet-induced potential and velocity. 
Substituting the approximation (B-75) into equation (B-3) and using hypothesis (B-6), 
we obtain: 




P 0 K1, 1)+M|I(2, 1) + ^ 1(1, 2)+^ 1(3, 1)+^ t? I(2,'2)+^m t?t? I(1,3). (B-85) 


Here 


I(M,N) = ^ 


r i (-* 

3(p 


n(M, N) -p(nq(M,N)-3p • pn(M,N)) 


(B-86) 
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where 


rf(M, N) = (-2aC(M + N), -2bC(M, N + 1),C(M,N)), 

nq(M, N) = -aC(M + 2, N) -bC(M, N + 2), 


(B-87) 

(B- 88 ) 


and 


pn(M, N) = (pT* ^(M + 1, N), p > -’n(M, N + 1), ap • rT(M + 2, N) + bp • n(M, N + 2)). (P-89) 

In computing the last component of pn the first two components* of 1? may ignored on 
account of hypothesis (B- 6 ). The induced velocity may be obtained by differentiating 
equation (B-85) and is given by 1 ° 

V = /U 0 T(1, 1) + M£?(2, l) + ju t7 T(l,2)+~^|T(3, 1 ) + ?(2, 2) + Jhrj V 1(1, 3)- O-90) 

Here 


J(M, N)^ 


-3 (n(M, N) - 3(p • n(M, N))$+ ((3nq(M, ff) r 1 5(p * pn(M, N))) p 


3pn(M, N) + 3np(M, N^) 


(B-91) 


where 


rip(M, N) = (-2ap • q(M + 1, N), ~2bp : cf(M, N f 1), p *~q(M, N^) . (B-92) 

In computing the first two components of np, the last component of q may be ignored on 
account of hypothesis (B- 6 ). Note that the quadrupole terip of expansion (B-75) is not 
used in the computation of doublet-induced potential and velocity. 

The computation of the C(M,N) integrals of equation (B-78) is similar to the 
computation of the H integrals of section (B.3) but much simpler. The range for the 
indices M and N is the same as that for the H integrals; i.e., 

M = 1, MXQ; N = 1,MXQ-M+1 (B-93) 

where MXQ is given in (B-36). Using the same potation as that of figure 31, we note 
that integration by parts yields 

,4 

C(M, N) ~ (M N ~- £ aG(M, N); M - 1, MXQ and N= l, MXQ - M + 1 (B-^4) 

where 

G(M, N) = f £ M_1 t ? N_1 dC; M = l,MXQand N= I.MXQ-M + 1 (Br95) 

J L 
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The integrals G(M,N) may be obtained with the aid of two recursion relations. We have 
the obvious identity • 

G(M + 1 , N) + v G(M, N + 1) = aG(M, N) (B-96) 

Integration by parts yields 

-(M-1)^G(M-1,N) + (N-1)^G(M,N-1) = D(M,N) (B-97) 

where 

D(M. N) = $ M_1 77 N_1 


' (B-98) 


The quantities D(M,N) may be evaluated directly or else recursively with the aid of 
equation (B-59). The recursion relations (B-96) and (B-97) may be recombined to yield 
an efficient procedure for evaluating the G integrals. - . \ * ■ 

Procedure 6: 


1. a. If li^l < \ug\ 


1 


(i) G(l, N) = D(l, N + 1); N = 1, MXQ. 

(ii) G(M, N) = -^2 g(M - 1,N+ 1) +— G(M - 1 , N); 
M = 2, MXQ and N = 1 , MXQ - M + 1 . 


(B-99) 


(B-100) 


b. If |^| < |^| 

. , , (i) G(M, 1) = D(M + 1,1); N = 1, MXQ 


v t a 

(ii) G(M, N) = — -G(M + 1, N- 1) + — G(M, N - 1); 


N = 2, MXQ and M = 1 , MXQ - N + 1 . 


(B-101) 


(B-102) 


This completes the evaluation of the source and doublet potential and velocity for a 
distant field point. 


87 



APPENDIX C 

GEOMETRY UPDATE COEFFICIENTS 


The coefficients dE/dO, dF/dd, and dGIdd of the Jacobian, equation (48), are germed 
geometry update coefficients. They are calculated based on the following assumptions: 


a) Changes of the aerodynamic influence coefficients due to changes pf the panel 
inclination 0 are neglected. 

b) The vector normal to the free sheet at control points close to wing edges is not 
affected by changes of 6. 

c) The chord length of panel segments in transverse geometry cuts (y,z-planps) dops 
not change during the iteration, i.e., li >m = constant (see fig. 32a). 


d) The free-sheet geometry is updated such that panel corner points move only in 
transverse cuts, i.e., Ax = 0. 


Coefficient dE IdO 


E = n * V is defined as the normal velocity at boundary points in the vicinity of wi 
edges. Hence, 




3E 

30 


3i? 

30 * 


V S + n • 


3 V s 
30 


= 0 


(C-l) 


The coefficient is zero due to assumptions a) and b) and because dn IdO = 0 at all 
boundary points on the wing. 


Coefficient dFIdO 


The quantity F stands for n • V° on the wing and the pres§ur$ juipp Acp ^cross the free 
sheet. Let 


F, = n • V" 


(wing) 


(C-2) 


and 


2 vS . 


F 0 = Acp =f — ^ V • • VM 
^ U ^ 

rir'i 


(free p^ieet) 

which is given by equation (E-ll). Application of assumption a) results in 


(C-3) 


3Fj _ ayS 
30 n 30 : 


(C-4) 
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3Fo 2 3 


(C-5) 


V/u. is a vector parallel to the free sheet; i.e., 


VM ’ n* = 0 


(C-6) 


Differentiating equation (C-6) for 6 gives 


9. , \ -*• _ _ 9if 

n - -VM- W 


(Cr7) 


Noting that d/50(V/u.)A0 is a vector normal to the free-sheet surface for small rotations 
A 0, and recalling that ii is the unit normal vector (x? • n = 1), equation (C-7) leads to 


^(VM) = -(w|g-)TT 


(C-8) 


Equation (C-5) thus takes the form 

9F-> 


90 


U„ 


2 l-f vS \/„ 9n 
i 11 ' V^JIVM 


(C-9) 


Expressions for dn/dd will be derived below. 

Coefficient dG/dO 

-►c 

G = n • V is the normal velocity component at the center of a free-sheet panel. Hence, 

9G = ^S.an (C-10) 

90 90 

taking assumption a) into account. 

Derivative dxt/dd 

The panel normal vector is 


n = 


N 


(C-ll) 


with 


illustrated in figure 32b. Hence, 


INI 

— ► — ► 

N = A x B 


9 n 
90 


INI 


9N U 9n\-> 
90 \ n ' 90 /" 


(C-12) 


9N _ 9A 9B 

90 gg x B + A x 90 


(C-^3) 
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The vectors A and B are expressed in terms of vectors pointing in direction of the panel 
corner points from some common origin. They read for panel i,m 


^i,m ^i,m + ^i+l,m ^i,m+l ^i+l,m+l 


(C-14) 


B 


i,m ^i,m + 1 ^i+ 1 ,m 1 ,m+ 1 


Assumptions c) and d) are used for the calculation of the changes in panel corner-point 
positions. Adopting the notation of^figure 32a for the mth transverse cut through the 
free sheet, the partial derivatives dP/d0 take the form 


3P 


i,m _ 


ae >,« 


3P 




90 J,n 


v — 


~ *j,n sin 0 ]n 
l j,n cos0 j,n 


(n ^ m or j i) 


(C-15) 


(n = m, j < i) 


Consequently, the derivatives of Aj >m , Bj >m , Nj >rn are zero for n^m, m + 1, or j 5= i. 
Hence, 


?I!ijn =Q 

M j,n 


(n ^ m, m + 1, or j > i) 
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APPENDIX D 

DOUBLET-STRENGTH UPDATE COEFFICIENTS 


The following coefficients of the Jacobian, equation (48), are termed doublet-strength 
update coefficients. 

3E 3E 3F 3F 3G 3G 
3m e ’ 1 3 Me 5 3 p r 5 3 m c ; 3 M r 


The boundary conditions F == 0 are divided into 


Fj = n • V s = 0 


Fo = Act 


(wing) 


= — — vS • 


V a • VM = 0 


tj — (free sheet/wake) 

^oo 

Coefficients 3E /d/x, dFi/d/x, dG/d/x (/x = /x e , fji T ) 

The symbols E, F, and G represent normal velocities at control points 

V f S = n • V s 

which can be written in subscripted notation as 

V ? S + u ? 

c 

The derivative of with respect to a particular doublet parameter /x n is 

dV ° 


jp— = Aj. 
op n 51 n 


(D-l) 


(D-2) 


(D-3) 


(D-4) 


(D-5) 


Hence, the coefficients dEldfj. e , dE/d/x r , dFj /d/j. e , 9Fi/3/x r , 3G/dp. e , dGldfi T are 
aerodynamic influence coefficients. 


Coefficients 3F2 /d/x (p. = /x e , /x r ) 

Differentiating F 2 with respect to a particular doublet parameter /x n gives 


9F- 




2 /dv b -»c a , \ 

VM + • — (Vju) 


Recalling equation (26) 


one obtains 


3p r 


V s = (A kj n, ♦ U k ) ? k (k = {,1,» 


av b 

3M n 


= A kn (k = T? > 0 


(D-6) 


(D-7) 


D-8) 
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Next, the gradient of — /x(£,t 7) is to be expressed in terms of the doublet parameters /xj 
before the indicated differentiation can be performed. Applying 

-j _ 9 — ► , 3 — ► , 3 _ -► 

V ‘ e r? 9f 

to the quadratic distribution of doublet strength defined by equation (10) yields 

VM = (M| + M^ + ^T?)e^ + (/i n + M^ + M w T?)e^ (D-9) 

The six coefficients of the panel doublet distribution are linearly related to the doublet 
parameters by 


M 0 


c oo 


r > 



C 0 £ 




► = 

1 

1 

1 

1 °° 

1 

1 

< 

"j 



c o^ 




C 0 ? 7 ? 



^nv 

J 


1 

n 

0 

-3 

L 


. J 


This is equation (41) rewritten with Coq, Cq^, Co^etc., denoting the rows of the matrix 
[Cq]. Then V^t takes the form 

V// = (c 0 £ + £C 0 ££ + vCq^ + (c Qr? + + J?C 07?T ^Mj} t v (D-ll) 

Hence, 

a^ ( VM) = ( C 0? + + ^ c 0^) n ^ + ( c 0r? + S C 0 £t? + 7 ? c 0r ? r?) n (D-12) 

Introducing the terms given by equations (D-7), (D-8), (D-ll), and (D-12) to equation 
(D-6) results in 
OF'} 9 / 

~dj^ = {[ A |n ( C 0£ + + vC m ) + A^ n (c 07? + + r?C 0r?7? )] j/ijf (D-13) 

+ v^(c 0| + sc m + nc 0 ^) n + v^ s (c 07? + |c 0 ^ + T?c 0 ^) n | 

The subscript n of 

A T?n ! ( C 0£ + «C(«€ + ^Oj-??)^ (^ 07 ? + ^0£t? + ^^07777^ n 
indicates the nth term of the row of the coefficients. 


93 



APPENDIX E 
PRESSURE EQUATIONS 


For incompressible flow the pressure coefficient cp reads 

-Hk — ► 

V • V 


Cp - 1 - 


(E-l) 


U, 


where V is the velocity vector and Uoo is the magnitude of the freestream velocity. The 
velocities of the two sides of a doublet sheet are different, say V^, on the upper side and 
Vj on the lower side. Introducing the average velocity V s and the jump in velocity 
V D , i.e., 




= t{Vv,) 


-►n — ► — ► 

v D = v u -v, 


one obtains 


V„ = ^S + ivD 
V, = \?S _^D 

Recalling that doublet strength /x is the jump in potential across the sheet, 

M = 4> u -^, 

the gradient of /u. 

Vju = V4> u - V<I>, = V u - Vj = V D 
is seen to be equivalent to the jump in velocity across the sheet. 


(E-2) 

(E-3) 

(E-4) 

(E-5) 

(E-6) 

(E-7) 


Hence, the pressure coefficients on upper and lower sides of the doublet sheet take the 
form 


C P, 


C P, 


1 - jj4(v S • \? S + V S • VM + JVM * Vvj 
1 — ■y ^ • Vm +^- Vm * Vt 


(E-8) 


(E-9) 
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The pressure jump, defined by 


becomes 



(E-10) 


Ac P = TT2 

u oo 


(E-ll) 
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